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the program is free from computational and logic errors, it cannot be considered validated. 
Any application of this program without additional verification is at the risk of the user. 



ABSTRACT 



This report presents a satellite simulation algorithm written for the IBM PC and clones 
in Microsoft QuickBASIC 4.5. The program can either generate an ephemeris of satellite 
ground-track position and location with respect to a fixed ground station, or determine the 
position and time the satellite is above the horizon of a specified ground station. Provision 
is made for the program to handle multiple satellites and multiple ground stations. Input 
is either by way of the computer keyboard or by separate files for satellites and ground 
stations. Output is written both on the screen and to an ASCII output file which can be 
accessed by a data based management program. 
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I. INTRODUCTION 



This report presents a satellite simulation algorithm written for the IBM PC and clones 
in Microsoft QuickBASIC 4.5. The program can either generate an ephemeris of satellite 
ground-track position and location with respect to a fixed ground station, or determine the 
position and time the satellite is above the horizon of a specified ground station. Provision 
is made for the program to handle multiple satellites and multiple ground stations. Input 
is either by way of the computer keyboard or by separate files for satellites and ground 
stations. Output is written both on the screen and to an ASCII output file which can be 
accessed by a data based management program. 

Two versions of the satellite kinematics subroutine pos. update are presented. These 
algorithms were provided by NAVSPASUR [Ref. 1], The first is a very accurate algorithm 
based upon a paper by Brouwer [Ref. 2] and is contained in NAVSPASURs FOR I RAN pro- 
gram called SHOWALL. A listing of the BASIC version of this algorithm is presented in 
Appendix E. The second algorithm is a much simplified approximation to the Brouwer 
model. Equations for the simplified model are contained in Appendix A, and the BASIC 
implemention is embedded in the listing of the full simulation program in Appendix D. 

The next sections contain program operating instructions, structure of the input/ 
output files, and suggestions for future improvements. The computation of the times of 
culmination, rise and set are in Appendix B, and the computation of satellite heading is in 
Appendix C. Readers unfamiliar with the nomenclature of satellite orbital elements should 
consult any standard work such as Escobal [Ref. 3, C'h. 3]. 

II. OPERATING INSTRUCTIONS 

A. Overview. The source code is named satsta.bas and the executable code is named 
satsta.exe. To run, enter SATSTA at the DOS prompt and press the <_j key. 

The first menu (Fig. 1) requests a data input mode. Pressing 1 or 2 will require the 
data for one satellite and one ground station to be input from the keyboard. Pressing 3 or 4 
assumes that two data files are accessible: either elements.dat or elements. nss containing 
the orbital elements for one or more satellites, and STATION.DAT containing the location of 
one or more ground stations — the structure of these files is discussed in the next section. 
No matter which of these four options is chosen, a simulation starting time, ending time 
and time-step increment will be required input from the keyboard in the third menu. 

The second menu (Fig. 2) requests a run mode. Option 1 creates an ephemeris of 
satellite geographic position and satellite position with respect to the ground station for 
equally spaced times for the entire simulation duration. Option 2 computes the satellite 
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geographic position and satellite position with respect to the ground station only for those 
times during the simulation duration at which the satellite is above the horizon of the 
ground station. 

The third menu (Fig. 3) requests the simulation-times, consisting of a starting time, 
an ending time, and a time-step increment (At). 

The fourth and fifth menus (Figs. 4 and 5, respectively) appear only if the first menu 
(data input mode) options 1 or 2 were chosen. Menu four requests the satellite epoch and 
elements, and menu five requests the station coordinates. 

No matter which options are chosen, output is directed to the screen and is appended 
to file output . DAT. If the file does not exist, it is created. This ASCII file is structured so that 
it can be used as input to a data base management program, if desired. The structure of 
this file is discussed in a later section. 

B. Input. Menu 1 requests the selection of a data input mode. Depending upon the 
source of satellite elements, there are two ways in which these data can be input. The 
first is using standard “epoch of date” elements in which the time of perigee passage is the 
epoch. The second is the NAVSPASUR ''One-Line Charlie” elements in which the longitude 
(or right ascension) of the ascending node has been referenced to 0000 hours, 1 January 
19S5. The ''One-Line Charlie” option has been added for those users who may obtain 
satellite elements from NAVSPASUR. An example of "epoch of date” input is shown in 
Fig. 3. 



Select data input mode 1, 2, 3 or 4: 

Keyboard input for satellite k ground station: 

1. Epoch of date. 

2. NAVSPASUR Reference Date. 

Data file input for satellite k ground station: 

3. Epoch of date. 

4. NAVSPASUR One-Line Charlie. 

Figure 1. Menu 1 — Select data input mode. 



If data file input Option 3 is chosen from Menu 1, two ASCII data files, elements.dat and 
station.dat, must have been previously created. The elements.dat file must contain a single 
line of for each satellite of interest. Each line must contain 10 entries, each enclosed in 
double quotations and separated by a blank space (Fig. 6). The entries are: (1) a satellite 
identification consisting of at most five alphanumeric characters; (2) the epoch date (date 
of perigee passage) in the form yyyy.mmdd, where yyyy is the year, mm is the month, and 
dd is the day; (3) the epoch time (time of perigee passage) in the form hh .mmss, where hh 
is the hour, mm is the minutes, and ss is the seconds; (4) the orbital eccentricity; (5) the 
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longitude (or right ascension) of the ascending node in degrees; (C) the orbital inclination 
in degrees: (7) the argument of perigee in degrees; (S) the mean anomaly in degrees: 
(9) the mean motion in revolutions per day; and (10) the orbital decay in revolutions per 
day-squared. 

If data file input Option 4 is chosen from Menu 1, two ASCII data files, elements. nss 
and station.dat, must have been previously created. The elements. nss file must contain a 
single line of for each satellite of interest in the NAVSPASUR format. Each line must contain 
10 entries starting in column 1 and ending in column 65 (Fig. 7). The entries are: (1) a 
satellite identification consisting of at most five alphanumeric characters (columns 1-5); 
(2) the mean anomaly in revolutions (columns 6-13); (3) the mean motion in radians per 
Herg 1 (columns 14-21); (4) the orbital decay in radians per Herg-squared (columns 22 
27); (5) the orbital eccentricity (columns 28-35); (6) the argument of perigee in revolutions 
(columns 36-43); (7) the longitude (or right ascension) of the ascending node in revolutions 
(columns 44-51): (S) the orbital inclination in revolutions (columns 52-59): and (9) the 
epoch date (date of perigee passage) in the form yymmdd, where yyyy is the year, mm is the 
month, and dd is the day (columns 60-65). A decimal point must be present in columns 6. 
14, 22, 2S. 36. 44 and 52. 



Select run mode 1 or 2: 

1. Create an ephemeris. 

2. Find times satellite is above the horizon. 

Figure 2. Menu 2 — Select run mode. 

The station.dat file must contain a single line of for each ground station or geographical 
location of interest. Each line must contain 5 entries, each enclosed in double quotations 
and separated by a blank space (Fig. S). The entries are: (1) a station identification 
consisting of at most three alphanumeric characters; (2) the station latitude in the form 
dd.mmss (minus if south), where dd is degrees, mm is minutes, and ss is seconds; (3) the sta- 
tion longitude in the form ddd.mmss (minus if west); (4) the station altitude in feet (minus 
if below sea level); and (5) the name of the station or station location (alphanumeric). 

If a data file input option is chosen, then SATSTA will open the two files described 
above and compute either an ephemeris or the times that each satellite is above the horizon 
(Option 1 or 2 of the second menu) for each ground station. Otherwise, keyboard input for 
one satellite and one ground station will be requested from the fourth and fifth menues, 
respectively. 

1 The Herg — a unit of time named after the astronomer Paul Herget — is defined as the period, divided by 
27 t, of a fictitious satellite at zero altitude over a smooth spherical earth with no atmosphere. 
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Simulation-Time Data: 



Starting date (yyyy.mmdd) 
Starting time (hh.mmss) 
Ending date (yyyy.mmdd) 
Ending time (hh.mmss) 

Time increment (minutes) 



? 1983.0201 
? 00.00 
? 1983.0202 
? 00.00 
? 01 



Press any key to continue. 

Figure 3. Menu 3 — Select simulation times. 



Five pieces of information regarding the simulation are requested from the third menu 
(Fig. 3): (1) the simulation starting date in the form yyyy .mmdd, (2) the simulation starting 
time in the form hh.mmss; (3) the simulation ending date in the form yyyy.mmdd; (4) the 
simulation ending time in the form hh.mmss; and (5) the simulation time-step increment 
in minutes. The screen will display sample data for each entry in the appropriate format. 
If any sample entry is to be changed, it should be completely overwritten in order for it 
to be read in properly. When each entry is correct, enter it into the computer by pressing 
the «_i key. 

If the keyboard entry option was chosen from the first menu, then two more menues 
will appear, otherwise the simulation will commence with input being obtained from the 
two data files. 



Satellite Data: 

Epoch date (yyyy.mmdd) ? 1983.0201 

Epoch time (hh.mmss) ? 00.00 

Input mean elements of epoch: 

Satellite ID Number ? 11111 

Eccentricity ? 0.0005545 

Long, ascending node (deg) ? 272.43497 

Inclination (deg) ? 65.06057 

Arg . of perigee (deg) ? 295.41470 

Mean anomaly (deg) ? 258.10682 

Mean motion (rev/day) ? 15.44194 

Decay (rev/day~2) ? 0 

Press any key to continue. 

Figure 4. Menu 4 — Input satellite data. 



Menu 4 (Fig. 4) requires ten separate entries: (1) the epoch date (date of perigee 
passage) in the form yyyy.mmdd, where yyyy is the year, mm is the month, and dd is the 
day; (2) the epoch time (time of perigee passage) in the form hh.mmss, where hh is the 
hour, mm is the minutes, and ss is the seconds; (3) a satellite identification consisting of at 
most five alphanumeric characters; (4) the orbital eccentricity; (5) the longitude (or right 
ascension) of the ascending node in degrees; (6) the orbital inclination in degrees; (7) the 
argument of perigee in degrees; (S) the mean anomaly in degrees; (9) the mean motion in 
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revolutions per day; and (10) the orbital decay in revolutions per day-squared. The screen 
will display sample data for each entry in the appropriate format. Again, if any sample 
entry is to be changed, it should be completely overwritten in order for it to be read in 
properly. When each entry is correct, enter it into the computer by pressing the _ key. 

Station Coordinates: 

Station Identification Number ? 001 

Station Location or Name ? Daisy, Tenn . 

Station latitude (dd.mmss, South minus) ? 35.12 

Station longitude (ddd.mmss, West minus) ? *85.12 

Station altitude (feet) ? 500.0 

Press any Xey to continue. 

Figure 5. Menu 5 — Input ground station data 

Menu 5 (Fig. 5) requires five separate entries: (1) a station identification consisting of 
at most three alphanumeric characters; (2) a station location or name: (3) the station lati- 
tude in the form dd.mmss (minus if south); (4) the station longitude in the form ddd.mmss 
(minus if west); and (5) the station altitude in feet (minus if below sea level). 

C. Output. Output is written to the screen and is appended to the data file output.dat 
if it exists, otherwise file output.dat is created. A sample from the output.dat file is shown 
in Fig. 9 — -the structure of this file is discussed later. 

Sample output from the “create an ephemeris" option screen is shown in Fig. 10. 
The apogee/perigee output is approximately valid only for the first orbit. In particular, 
the longitude is not corrected for the earth’s rotation during the current orbit. Also, the 
ascending node and argument of perigee axe not corrected for their rates of change during 
the current orbit. It should be further noted that for nearly circular orbits, the height of 
the satellite will probably get lower than the perigee value and higher than the apogee 
value at locations other than the perigee and apogee positions — this is an effect of the 
earth’s oblateness. 

The main body of the ephemeris output consists of date, time, satellite latitude, 
longitude and height, followed by the elevation angle, azimuth, and range (in kilometers) 
of the satellite from the station. The satellite “look angle' 1 is the angle between the satellite 
“ground-zero” point (usually called the satellite “sub-point”) and the station. The last 
column of output is the satellite heading. 

A sample of the output from the “times satellite is above the horizon" option screen 
for a time-step (At) of 1 minute is shown in Fig. IF First the time of culmination for 
the current orbit is computed and saved. The culmination time is then incremented and 
decremented by At until the satellite is below the horizon. The time of rise and set is 
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Figure 9. Sample output.dat file (some output expunged). Epoch of date. 



Satellite: 11111 



Station: 001 



Daisy, Tenn. 



Apogee : 
Height = 

Latitude = 
Longitude = 

Perigee : 
Height = 

Latitude = 
Longitude = 



442.8 km. 

-55.1540 

100.2099 



450.4 km. 

55.1539 

280.2099 



Tue 1 Feb 1983 



year 


mo 


da 


hr 


mn 


sec lat 


long 


height 

(km) 


elev 


azim 


range 


look 

angle 


head 


1983 


2 


1 


0 


0 


0.0 -12.26 


327.56 


434.0 


-31.60 


123.30 


7437 


52.90 


157.82 


1983 


2 


1 


0 


10 


0.0 -45.84 


347.68 


441.5 


-50.59 


136.39 


10382 


36.56 


145.02 


1983 


2 


1 


0 


20 


0.0 -65.19 


48.39 


446.4 


-69.07 


151.85 


12354 


19.68 


83 .46 


1983 


2 


1 


0 


30 


0.0 -44.71 


106.99 


439.1 


-83.25 


221.39 


13089 


6.45 


34.07 


1983 


2 


1 


0 


40 


0.0 -10.95 


126.55 


430.7 


-70.51 


302.57 


12481 


18.11 


22.04 


1983 


2 


1 


0 


50 


0.0 24.08 


141.07 


435.9 


-51.76 


316.67 


10574 


35.26 


24.38 


1983 


2 


1 


1 


0 


0.0 55.55 


168.97 


449.2 


-32.06 


324.49 


7561 


52.31 


46.47 








F 


igure 10. Output from ephemeris option. 


Epoch of date. 







Satellite : 


11111 


Station 


: 001 


Daisy, Tenn. 


























height 








look 




year 


mo 


da 


hr 


mn 


sec 


lat 


long 


(km) 


elev 


azim 


range 


angle 


head 


1983 


2 


1 


1 


13 


19.4 


54.37 


263.10 


450.1 


0.00 


340.48 


2436 


69.23 


135.43 


1983 


2 


1 


1 


13 


42.1 


53.35 


264.76 


449.7 


1.39 


341.68 


2285 


69.19 


136.92 


1983 


2 


1 


1 


14 


42.1 


50.53 


268.75 


448.7 


5.55 


345.74 


1892 


68.57 


140.46 


1983 


2 


1 


1 


15 


42.1 


47.58 


272.27 


447.6 


10.78 


351.87 


1514 


66.79 


143.50 


1983 


2 


1 


1 


16 


42.1 


44.51 


275.37 


446.4 


17.70 


2.15 


1167 


63.08 


146.09 


1983 


2 


1 


1 


17 


42.1 


41.35 


278.14 


445.1 


26.48 


21.39 


892 


56.95 


143.32 


1983 


2 


1 


1 


18 


42.1 


38.12 


280.62 


443.8 


32.24 


55.33 


773 


52.35 


150.23 


1983 


2 


1 


1 


19 


42. 1 


34.83 


282.88 


442.6 


27.03 


90.21 


875 


56.40 


151.86 


1983 


2 


1 


1 


20 


42.1 


31.49 


284.94 


441.3 


18.10 


110.45 


1141 


62.69 


153.26 


1983 


2 


1 


1 


21 


42.1 


28.11 


286.84 


440.1 


10.97 


121.17 


1484 


66.60 


154.46 


1983 


2 


1 


1 


22 


42.1 


24.69 


288.61 


438.9 


5.61 


127.48 


1862 


68.51 


155.46 


1983 


2 


1 


1 


23 


42.1 


21.25 


290.27 


437.9 


1.36 


131.58 


2255 


69.20 


156.31 


1983 


2 


1 


1 


24 


3.8 


19.99 


290.85 


437.5 


-0.00 


132.74 


2400 


69.25 


156.58 



Figure 11. Output from ‘‘times above horizon” option. A / = 1 minute. Epoch of date. 



iterated until the time at which the satellite is on the horizon. All of these times and 
positions are then output to the screen. Note that the information printed out in this 
option, with the exception of the apogee/perigee positions, is identical to that printed 
out in the “ephemeris” option. By making At large, say 10 minutes, in the kC times above 
horizon'’ option, then only the times of satellite rise, culmination, and set are computed 
and output (see Fig. 12). 

Each line in the output.dat file consists of the station identifier and the satellite iden- 
tifier followed by the same data, and in the same order, as that printed to the screen. This 
can be seen by comparing the first six lines of Fig. 12 with the same lines in Fig. 9. 



i 



Satellite: 11111 Station: 001 Daisy, Tenn. 



height look 



year 


mo 


da 


hr 


mn 


sec 


lat 


long 


(km) 


elev 


azim 


range 


angle 


head 


1983 


2 


1 


1 


13 


19.4 


54.37 


263.10 


450.1 


0.00 


340.48 


2436 


69.23 


135.43 


1983 


2 


1 


1 


18 


42.1 


38.12 


280.62 


443.8 


32.24 


55.33 


773 


52.35 


150.23 


1983 


2 


1 


1 


24 


3.8 


19.99 


290.85 


437.5 


-0.00 


132.74 


2400 


69.25 


156.58 


1983 


2 


1 


2 


49 


40.9 


45.38 


250.95 


446.7 


0.00 


306.66 


2428 


69.23 


145.41 


1983 


2 


1 


2 


54 


12.9 


30.63 


261.84 


441.0 


11.96 


251.45 


1429 


66.15 


153.58 


1983 


2 


1 


2 


58 


44.7 


15.05 


269.44 


436.2 


-0.00 


195.12 


2394 


69.26 


157.45 


1983 


2 


1 


14 


47 


13.6 


16.76 


285.33 


433.6 


0.00 


150.29 


2387 


69.32 


22.83 


1983 


2 


1 


14 


50 


58.0 


29.64 


291.62 


437.9 


6.44 


106.44 


1792 


68.35 


26.06 


1983 


2 


1 


14 


54 


43.3 


42.02 


300.04 


443.2 


-0.00 


63.10 


2419 


69.28 


32.12 


1983 


2 


1 


16 


20 


48.3 


17.88 


262.22 


433.9 


0.00 


216.51 


2389 


69.32 


23.02 


1983 


2 


1 


16 


26 


16.7 


36.50 


272.26 


440.8 


55.22 


303.76 


528 


32.33 


28.94 


1983 


2 


1 


16 


31 


46.6 


53.37 


289.30 


448.2 


-0.00 


24.72 


2431 


69.26 


43.11 


1983 


2 


1 


17 


59 


48.3 


37.38 


249.28 


441.2 


0.00 


283.77 


2414 


69.29 


29.38 


1983 


2 


1 


18 


2 


51.4 


47.05 


257.60 


445.5 


3.50 


317.85 


2066 


69.01 


36.02 


1983 


2 


1 


18 


5 


55.6 


55.73 


269.75 


449.2 


-0.00 


351.87 


2433 


69.26 


46.77 



Figure 12. Output from “times above horizon' 1 option. A t = 10 minutes. Epoch of date. 



Fig. 13 is a sample output obtained using the NAVSPASUR reference time. A compar- 
ison of Fig. 10 and Fig. 13 shows that the satellite latitude, height, and heading is the 
same but the longitude differs because of the difference of the earth’s rotation between the 
two different reference times. The satellite-station relationships differ for the same rea- 
son. Figures 12 and 14 are completely dissimilar because of the differing satellite-station 
relationships. 



Satellite: 11111 Station: 001 Daisy, Tenn. 



Apogee : 
Height = 

Latitude = 
Longitude = 

Perigee : 
Height = 

Latitude = 
Longitude = 



442.8 km. 

-55.1540 

230.8483 



450.4 km. 
55.1539 
50.8483 



year 


mo 


da 


hr 


mn 


sec 


lat 


long 


height 

(km) 


elev 


azim 


range 


look 

angle 


head 


1983 


2 


1 


0 


0 


0.0 


-12.26 


98.20 


434.0 


-77.96 


351.78 


12919 


11.09 


157.82 


1983 


2 


1 


0 


10 


0.0 


-45.84 


118.32 


441.5 


-79.40 


232.53 


12965 


10.02 


145.02 


1983 


2 


1 


0 


20 


0.0 


-65.19 


179.03 


446.4 


-60.98 


210.44 


11620 


27.15 


88.46 


1983 


2 


1 


0 


30 


0.0 


-44.71 


237.62 


439.1 


-41.32 


205.69 


9023 


44.80 


34.07 


1983 


2 


1 


0 


40 


0.0 


-10.95 


257.18 


430.7 


-20.38 


203.40 


5458 


61 .42 


22.04 


1983 


2 


1 


0 


50 


0.0 


24.08 


271.70 


435.9 


12.73 


194.92 


1373 


65.77 


24.38 


1983 


2 


1 


1 


0 


0.0 


55.55 


299.61 


449.2 


-5.07 


32.07 


3062 


68.64 


46.47 



Figure 13. Output from ephemeris option. NAVSPASUR reference time. 
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Satellite : 



11111 



Station: 001 Daisy, Tenn. 



year 


mo 


da 


hr 


mn 


sec 


lat 


long 


height 

(km) 


elev 


azira 


range 


look 

angle 


head 


1983 


2 


1 


0 


47 


32.1 


15.54 


267.79 


433.2 


0.00 


199.85 


2385 


69.32 


22.63 


1983 


2 


1 


0 


52 


59.8 


34.24 


277.36 


439.8 


58.70 


112.56 


509 


29.01 


27.87 


1983 


2 


1 


0 


58 


28.6 


51.41 


293.03 


447.4 


-0.00 


33.16 


2429 


69.27 


40.56 


1983 


2 


1 


2 


24 


57.1 


29.86 


251.15 


438.0 


0.00 


262.00 


2405 


69.30 


26.14 


1983 


2 


1 


2 


29 


10.7 


43.72 


260.91 


443.9 


8.69 


312.74 


1641 


67.65 


33.31 


1983 


2 


1 


2 


33 


25.7 


55.99 


276.83 


449.3 


-0.00 


3.07 


2433 


69.26 


47.22 


1983 


2 


1 


8 


58 


59.6 


56.01 


272.80 


450.7 


0.00 


356.74 


2437 


69.23 


132.74 


1983 


2 


1 


9 


3 


16.0 


43.68 


288.80 


446.0 


8.83 


47.32 


1637 


67.55 


146.72 


1983 


2 


1 


9 


7 


31.7 


29.71 


298.62 


440.7 


-0.00 


98.40 


2412 


69.24 


153.92 


1983 


2 


1 


10 


33 


57.2 


51.41 


256.65 


449.0 


0.00 


326.74 


2433 


69.23 


139.44 


1983 


2 


1 


10 


39 


26.8 


34.21 


272.34 


442.3 


58.17 


247.39 


514 


29.48 


152.14 


1983 


2 


1 


10 


44 


55.7 


15.45 


281.93 


436.3 


-0.00 


160.33 


2394 


69.26 


157.39 


1983 


2 


1 


12 


13 


31.5 


31.47 


250.42 


441.3 


0.00 


266.88 


2414 


69.24 


153.27 


1983 


2 


1 


12 


14 


27.8 


28.30 


252.20 


440.2 


0.29 


256.89 


2378 


69.24 


154.39 


1983 


2 


1 


12 


15 


23.2 


25.15 


253.85 


439.1 


-0.00 


247.03 


2407 


69.25 


155.34 



Figure 14. Output from “times above horizon" option. — 10 minutes. Navspasl'R reference time. 



III. SOURCES 

As mentioned in the introduction, the high resolution and simplified satellite kinematics 
models were obtained from NAYSPASUR [Ref. lj. The high resolution model is based upon 
a paper by Brouwer [Ref. 2], The code for the high resolution model is given Appendix E. 
while the code for the simplified model is contained in the pos. update subroutine in 
Appendix D. Both of these routines require a solution of Kepler's equation. A non-iterative 
state-of-the-art solution of Kepler’s equation, based upon a paper by Mikkola [Ref. G] is 
contained in subroutine kepler listed in Appendix D. 

The ground track determination is essentially that given by Escobal [Ref. 3. C’h. 10.3.2] 
but with a minor improvement described by Shudde [Ref. 7j. 

The method of computation of the times of satellite culmination, rise, and set are given 
in Appt lix B, and the computation of the satellite heading is contained in Appendix C. 
Additional references are given there. 

The IAU (1976) System of Astronomical Constants [Ref. 10, pg. S7] is used in SATSTA 
as well as the J2000.0 reference frame [Ref. 10, pg. Ll]. Satellite positions computed by 
NAVSPASUR’s SHOWALL [Ref. l] use the WGS 72 constants and the J1900.0 reference frame. 
For this reason, the agreement between SATSTA and SHOWALL is not exact. 



IV. DISCUSSION 

Although the position updating in SATSTA is based upon the satellite kinematics model 
of Brouwer [Ref. 2], other models are currently in use. Most notable is the work of Hoots 
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and Roehricli [Ref. 5]. SATRAK [Ref. 4] is a FORTRAN program implementing their models 
and is available from the Aerospace Defense Command, Peterson AFB. Unfortunately, the 
source code is not available for modification. 

The program presented here can be used as a basis for a more sophisticated model. 
The input scheme could be improved by using an input editor to catch improper types of 
entry and to allow corrections or changes to be made to previous entries. File management 
programs can be used to select satellite elements and station parameters from master files 
and insert them into “working” files for program input. 

Various sensor types have not been modelled and no provision has been made for 
including them in the program. Also, satellite-satellite visablity has not been modelled. If 
such features are desirable, they could be included by modifying the existing program. 
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APPENDIX A. THE SIMPLIFIED KINEMATICS ALGORITHM. 



The computation of the perturbed satellite orbit about the earth is performed using a 
method which is an approximation of that used by NAVSPASUR [Ref. 1], According to 
NAVSPASUR, this approximation will yield satellite positons with a random error of about 
10-15 nautical miles for a 20 day interval around epoch. 

The mean elements at epoch, which must be input, are: 

Eccentricity eo 

Right ascension of the ascending node Qq 

Inclination io 

Argument of perigee u 0 

Mean anomaly A/ 0 

Mean motion A/j , and 

Decay A /2 

Using the mean elements, the following auxiliary quantities are computed: 



- if" 2 / 3 



a 0 = A/j 



a 0 — a 0 



1 + 



C(3 cos 2 io — 1 ) 



2/3 



Got 1 ~ e l) 



2 rf /2 



Co — c 



n = -2 c 



5 cos' 1 io — 1 

« 2(1 - 4? ’ 

cos i 0 



a — — 



°o(l — e o) 2 
4a 0 



and 



3A/i 



A/ 2 , 



(AM) 



(A la) 

(A -2) 
(A 3) 
(A -4) 



where C = 3 /t J 2 is related to the oblateness of the earth. NAVSPASUR states the following: 

For optimal accuracy, ao from Equation A-l is used in Equation A-la to obtain a final 
value for a q. This latter ao is then used throughout. Where program size is a major 
consideration and additional errors of some ten miles in satellite height can be tolerated 
the user may choose to leave out Equation A-la and/or Equation A-4. 

Using the mean values of the elements at epoch To and the results of the above 
calculations, the mean values of the elements for any time t axe: 



M = Mo + M x (t - To) + M 2 (t -To) 2 . 

€ = Cq , 

uj — cjq + — A/q)) 
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£1 — Q 0 + (l(M ~~ Mo) ~ u e (t — To), 

% “ %o i and 
a = a 0 + a(t - T 0 ), 



where u> e is the earth’s sidereal rotation rate. 
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APPENDIX B. THE TIME OF CULMINATION, RISE AND SET ALGO- 
RITHM. 



In this program the earth is considered to be non-rotating, i.e., the coordinates of a station 
are given by latitude and longitude with the longitude unaffected by the earth’s rotation. 
Rotation is accomplished by subtracting the accrued rotation from the right ascension of 
the ascending node of the satellite. Although this departs from realism, the benefit is that 
the rectangular coordinates of the station need to be calculated only once. 





Fig. B-l. Satellite/Station Geometry. 



Referring to Fig. B 1, culmination of the satellite, S, with respect to a ground station, 
P, occurs when the elevation angle, /?, of the satellite is a maximum, which is equivalent 
to the minimization of the of the zenith angle, 90° — /z, which is also equivalent to the 
minimization of the angle // between z, the unit vector normal to the station location, and 
r, the vector from the earth’s center, 0, to the satellite. The equation for 7] is 



z r 

COS 71 = 

r 



(B-l) 



where z is given by 



z = 



cos (f) cos A 
cos (p sin A , 
sin (f> 



( B— 2 ) 



and 4> and A are the station geodetic latitude and longitude, respectively. The satellite 
position vector r is obtained from [Ref. 3, Equ. (3.57)] 



r = x^P + y w .Q, 



(B-3) 
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where x u and y u measure the satellite position in the orbital plane. 1 P and Q are the 
first two columns of the Euler matrix [Ref. 3, Equs. (3.43) and (3.44)] — they transform the 
in-plane x u and y ^ coordinates to the equatorial rectangular system and are given by 





'Px 




cos uj cos r — sin uj sin T cos i 




p = 


Py 


= 


cos LJ sin T + sin uj cos V cos i 


(B-4a) 




P: 




sin uj sin i 





and 





'Qx' 




— sin uj cos r — cos uj sin T cos i 


Q = 


Qy 


= 


— sin uj sin T + cos uj cos T cos i 




Qx 




cos uj sin i 



(B-4b) 



where u; is the argument of perigee and i is the orbital inclination; T = ft — u e (t — f 0 ), 
where ft is the longitude of the ascending node at the time, to, of the last element update; 
u e is the earth’s sidereal rotation rate and t is the current time. 

The components of r vary with time. We will transform the equations in such a 
way that will replace time by the eccentric anomaly E as the independent variable. The 
inclination i is time independent and the elements ft and u> will be considered to be 
independent of time over the period of a single orbit even though they are time dependent 
through perturbations. The components of r which will be considered to be time dependent 
are then i w , y w and T. 

Write Equation (B-l) as 



z • r — r cos rj = 0. 



Take the derivative with respect to E to obtain, 



z • 



dr 

dE 



dr 



drj 



— cos, + rs, nI) — 



0 . 



To find the value of E which minimizes rj, set drj/dE = 0. Thus, 



dr 

dE 



dr 

— cos , = °, 



or 




1 dr 
rdE Z 



r = 0. 



Since this equation cannot be solved explicitly for E, set 



F(E) = z 



dr 

dE 



1 dr 
rdE 



r 



(B-5) 



lr The origin of the system is at the focus of the orbit. X ^ is positive in the direction of the perifocus. y ^ 
lies in the orbital plane and is orthogonal to 
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and find E such that F(E) = 0 by the Newton- Raphson method. The Newton- Raphson 
procedure also requires E'(E), the derivative of F(E) with respect to E: 



F'(E ) = z 



d 2 r 1 / dr \ 
dE 2 + F \dE/ 



1 d 2 r 1 dr dr 

z • r — z • . 

r dE 2 r dE dE 



The first and second derivatives of Equation (B-3) are 



dE dE P + XuJ dE + dE®^ V ~~ dE 



and 

d 2 r _ d 2 x^ ^dx^dP d 2 P 0 dy w dQ d 2 Q 

dE 2 dE 2 + “ dE dE + ^dE 2 + dE 2 ^ + “ dE dE + ^ dE 2 ' 



(B-C) 



(B-7a) 



(B-7b) 



where [Ref. 3, Equ. (3.G9)] 

x u — a ( cos E — e) and 
y^ = a\/l — e 2 sin E. 

with first and second derivatives 



dx w 

7e 



—a sin E. 



dy^. 

dE 




e 2 cos E. 



and 



d 2 y^ 

dE 2 

d 2 Xj 

dE 2 



— a\/l — e 2 sin E = — y^, 
—a cos E. 



The first and second derivatives of Equations (B-4) are 



dP z 

dE 

dPy 

dE 

dP z 

dE 

dQ z 

dE 

dQy 

IF 

dQ. 

dE 



— Fy , ■ 



d£ 

'dE' 



= P d S 

1 dE' 

= 0, 

d£ 

“ ~ Qy dE' 

= Q - 

Wl dE' 



= 0 , 



(B-Saj 



(B-Sb) 



( B -Sc J 



(B-9a) 
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and 



2 






dE 2 

d 2 P, 



dE 2 

d 2 P z 



i._ a ,4£V + p. 



y \dE J 



d 2 r 

dE 2 ' 

d 2 r 

d£ 2 ’ 



d£ 2 

d 2 Q 



d£ 2 



= 0, 

= -Q, 



x ^ , dry d 2 r 

- 1 ^£2 



(B-9b) 



d£y 



d 2 r 



^ = _n (*£) +Q 
d £ 2 + Vl d £ 2 

d 2 Q z 

dE 2 



and 



= 0. 



Also required is [Ref. 3, Equ. (3.67)] 



and it’s derivatives 



r = a(l — e cos E) 
dr 



dE 

d 2 r 

dE 2 



= ae sin £, and 
— ae cos E. 



(B-lOa) 



(B-lOb) 



Finally, using Kepler's equation [Ref. 3. Equ. (3.79)], 

E — e sin E 



t = 



n 



+ T, 



(B-ll) 



where n is the rate of orbital motion and T is the time of last perifocal passage, the 
equation for T can be written as 



„ ~ , E — e sin E m . 

r = £1 — oj e [ (- T — to ) • 



n 



Then, 



and 



dr 

dE 



= -w. 



1 — e cos E ' 



n 



(B-12a) 



(B-12b) 



d 2 T — u; e (e sin E) 



dE 2 n ' ^ 12C ^ 

To solve Equation (B-5), an initial estimate of E is made and the right-hand side 
of Equation (B-13) is evaluated. The result, the left-hand side of Equation (B-13), is an 

improved estimate of £, 

F(E) 



E = E 



F'(EY 



(B-13) 
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This iteration is continued until \F(E)/ F'(E)\ is less than some prescribed amount. The 
value of E obtained will either maximize or minimize 77. The desired value of 77 will be the 
one for which cos 77 > 0. 

Once the eccentric anomaly for culmination, E c , is found, the time of culmination is 
determined directly from Equation (B-ll). Perturbing E c by plus or minus about 10°, it 
is possible to obtain good starting values for the determination of the satellite set and rise 
times, respectively. 

The procedure for finding the rise and set times is very similar to finding the time 
of culmination. The rise/set equations are developed by Escobal [Ref. 3, Sec. 5.4], but 
the algorithm presented here is somewhat simplified since, unlike Escobal, the time of 
culmination is available. In fact, many of the same equations can be used. Again, the 
eccentric anomaly will be used as the independent variable. Referring to Fig. B-l, satellite 
rise or set occurs when the elevation angle h is equal to zero, where h can be found from 
the equation 

sin h —— — , (B-14) 

P 



where 



p = r — R. 



( B— 15) 



In Equation (B-15), R is the station coordinate vector and p is the slant range vector 
from the station to the satellite. For the spheroid earth [Ref. 3, Equ. (5.54)]. 



R = 



G 1 cos (pcos A 
G 1 cos <f> sin A 
G 2 sin <f> 



(B-1G) 



where <t> and A are the station geodetic latitude and longitude, respectively. G\ and G 2 
are factors which account for the earth’s oblateness, and are given by 



G\ = — = + H and 

Y 1 — (2/ — / 2 ) sin 2 <j) 

= V-f 2 )?* + H, 

(2/ ~ P) sin 2 



(B-l 7 ) 



where 



a e = earth’s equatorial radius, 

/ = earth’s flattening factor, and 
H = station elevation above the ellipsoid. 
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Write Equation (B-14) as 

p ■ z = p sin h. 

When the satellite is on the horizon, h = 0 and sin h = 0. Since the resulting equation, 
p • z = 0, cannot be solved explicitly for E, set 

G{E)=p- z 



and find E such that G(E) = 0 by the Newton-Rap hson method. Using Equation (B-15), 
write G(E) as 

G(E) = z • (r — R). (B-18) 

The Newton-Raphson procedure also requires G'(E), the derivative of G(E) with respect 
to E. Since the earth is being treated as non-rotating with the sidereal rotation subtracted 
from the satellite’s ascending node. R is independent of E , hence 



G'(E) 




(B-19) 



In Equations (B-18) and (B-19), z, r, dr/dE, and the additional needed derivatives have 
all been given previously in this section. 

To solve Equation (B-18), and initial estimate of E is made and the right-hand side 
of Equation (B-20) is evaluated. The result, the left-hand side of Equation (B-20), is an 
improved estimate of E , 



E = E - 



C(E) 
G'(E )' 



(B-20) 



The iteration is continued until \G{E) / G' {E)\ is less than some prescribed amount. The 
corresponding time of rise or set is then obtained from Equation (B-ll). 
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APPENDIX C. THE SATELLITE HEADING COMPUTATION. 

A satellite can be located by a position vector x = xii -f x 2 j + rak in a rectangular 
coordinate system with origin at the earth’s center. In matrix form, and with Xj, x 2 and 
X3 expressed in spherical coordinates, 



x = 



r cos d> cos A 
r cos <p sin A 
r sin 



where <f> is the latitude, A is the longitude at the point, and r = |x| is the distance of the 
point from the earth’s center. 

Let z be the unit vector in the direction of r. Then 



z = 



-1 




cos 0 cos A 


z 2 


== 


cos q sin A 


. Z 3 _ 




sin <p 



It is convenient to define two other unit vectors in the plane tangent to z. These are local 
north n and local east e, defined by 



n = 



M 



— sin 0 cos A 






— sin A 


— sin 0 sin A 


and 


za 

e = . . = 


cos A 


COS <j) 




za 


0 



where 



- dz 1 

z^- — and 



dz 

Za = d\ ' 



d<p 

The vectors z, n and e provide the basis for a right-handed orthogonal coordinate system. 

If a satellite’s velocity vector v = tq i -f e 2 j + V3U. is known at position x then it is 
easily shown [Ref. 8] that 

n • v = cos a and 
e • v = sin a, 

where a is the course angle or heading. From these relations, a is found to be 

a = qatn(e • v, n • v) 

where qatn(y,x) is arctan(y/x) corrected to the proper quadrant. 

Since the position and velocity component’s are in rectangular coordinates, it will be 
necessary to convert them to spherical coordinates in order to determine the components 
of n and e. However, the spherical coordinate step can be circumvented as follows: First, 



cos q = n • v = — t'i sin <p cos A — t> 2 sin o sin A + t’3 cos C>. 
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Multiply both sides by cos <f> to give 



cos <j) cos a = — t)j sin <f> cos <j> cos A — v 2 sin <j> cos <j> sin A + u 3 cos 2 <f> 

= — sin <t>(z\vi + Z2V2) + V3 cos 2 <f> + Z3V3 sin <}> — Z3V3 sin <f> 
= - s\n(j>(ziVi + Z2V2 + Z3V3) + v 3 cos 2 <f> -f Z3V3 sin<^» 

= —Z3(z • v) + V3 cos 2 <f> + v 3 sin 2 <j> 

= V3- Z 3 ( Z • V). 



Siinilaxly, 



so that 



Then, 



sin a = e • v = — Vi sin A + i >2 cos A, 

cos (*> sin a = V 2 cos <j> cos A — Vj cos <p sin A 
= Zl v 2 - z 2 v x . 



a = qatn(cos <f> sin a, cos <j> cos a). 



Since cos <f> > 0 for all <f> 6 ( — 7r/2, 7r/2), the cos<^> term effectively ‘'cancels out” [Ref. 9] 
and so the heading is determined by 



a = qatn(sin a, cos a) 

= qatn(r 1 v 2 - -2*9, - z 3 (z ■ v)). 

This equation is used by NAVSPASUR [Ref. 1] in SHOWALL. In SATSTA, the equivalent form, 



a = qatn 



( 



X1V2 — x 2 vi 



V3 - 



X 3 (x • v) 



is used. 
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APPENDIX D. THE SATSTA PROGRAM LISTING. 



'SATSTA 

' Satellite Track. 03-03-88. Rev. 09-26-89 0 0830. 



DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 



FUNCTION acos# (xt) 

FUNCTION asin# (xf) 

FUNCTION at an2# (yf, xf) 

FUNCTION decimal# (v$) 

FUNCTION dmod# (xt, mf) 

FUNCTION dot# (a#(), b#()) 

FUNCTION gmst# (tu#, tm#) 

SUB apo. peri. gee (elem#(), elemufQ) 

SUB culmnate (zv#(), ve#, tO#, eat, coseta#, elemu#()) 

SUB euler (peri#, node#, incl#, psn#(), vel#()) 

SUB ground. track (r#, sat.dc#, sat. It#, sat.ht#) 

SUB hr. min. sec ( t ime . of . day# , hr % , minX, sec#, nX) 

SUB jd. to. date (jd#(), mk , d* , ylt, h#, mo$, day$) 

SUB j ulian . day . number (m&, d& , y k, ut#, jd#(), dj#, tj#) 

SUB kepler (m# , ec#, eat) 

SUB output4 (k*/,, yr&, mo&, da&, t ime . of . day# , It#, In#, ht#, el#, 
az#, rng#, langle#, head#, sta.idS, sat.idS) 

SUB pos. update (1*/., elem#(), time#, sat.x#(), sat.xd#(), elemu#(), 
r# , ea#) 

SUB rise. set (zv#(), sta.R#(), we#, tO#, ea# , elemu#()) 

SUB sat. head (r#, sat.x#(), sat.xd#(), head#) 

SUB sat . sphere . coord (sat.x#(), r# , tj#, tm#, sat.dc#, sat. In#, sat.ra#) 
SUB sta.p. coord (It#, In#, ht#, x#()) 

SUB sta.sat (sta.x#(), sat.xtQ, rng#, az#, el#, langle#) 

SUB time. update (k*/., deltat#, tm# , jd#, jd#(), tj#, t ime . of . day# , 
m k, dk , yk, h# , mo$, day$) 

SUB yrmoday (v$, yr&, mo&, dy&) 



* COLOR 7 , 1 
' SCREEN 0, 1 



DEFDBL A-Z 
OPTION BASE 1 

DIM jdepoch(2), jdstart(2), jdend(2), jd(2), sat.x(3), sat.xd(3), sta.x(3, 4) 
DIM elem(9), elemu(9), pv(3), qv(3) , sta.z(3), sta.R(3), jd . charli e(2) 

1 constants 

CONST pi = 3.141592653589793#, rad = pi / 180#, tvopi = pi + pi 
CONST j 2 = .0010826318# 

CONST nmi.ft = .3048# / 1852# ’ n.rai/foot 

CONST ae = 6378.14# ’ earth's equatorial radius (km.) 

CONST GE = 398600.5# ' g©0- grav . const. (km)*3/(sec)"2 

k = 60# * SQR(GE / ae) / ae ' ( eru) ~ (3/2)/min 
Herg.min = k 

CONST fl = 1# / 298.257# ' earth’s flattening factor 

CONST e.sqr = (2# - fl) * fl ' earth's ellipticity squared 

CONST we = 1.002737909350795# * tvopi / 1440# ' earth ' 8 rotation (radians/minute) 
CONST sixty = 60# 

' Sample input (NORAD Project Spacetrack Rpt . #3, Dec. 1980): 
amajorS = "1.04050189": ecc$ = "0.0086731": aincl$ = "72.8435" 
per i$ = "52.6988": anode$ = "115.9689": manomaly$ = "110.5714" 
motions = "16.05824518": decayS = "0" 
epochdateS = "1980.1001": epochtimeS = "23.41241138" 
startdateS = "1980.1001": starttimeS = "23.41241138" 
enddateS = "1980.1003": endtimeS = "0000" 
deltatS = "4" 



' Sample input of arbitrary higher altitude, higher ecc orbit 
ecc$ = "0.0200000": ainclS = "64.94868" 

periS = "202.277268": anodeS = "55.552464": manomalyS = "126.670896" 
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motion* = "12.00000": decay* = "4 . 281098d-3 " 
epochdate* = "1983.0201": epochtime* = "00.00" 
startdate* = "1983.0201": starttime* = "00.00" 
enddate* = "1983.0208": endtime* = "00.00" 
deltat* = "15" 

sta.lat* = "35.12": sta.lon* = "-85. 12": sta.alt* = "500.0" 

1 Sample input (NAVSPASUR: Sat. 22222 data from SH0WALL, private comm. 

* Mr. Fred Lipp, July 1988): 
sat. id* = "22222" 

ecc* = "0.0023064": aincl* = "64.94868" 

peri* = "202.277268": anode* = "55.552464": manomaly* = "126.670896" 

motion* = "16.06302": decay* = "4 . 281 1 19d-3" 

epochdate* = "1983.0201": epochtime* = "00.00" 

startdate* = "1983.0201": starttime* = "00.00" 

enddate* = "1983.0203": endtime* = "00.00" 

deltat* = "01" 

sta.id* = "001": sta.name* = "Daisy, Tenn." 

sta.lat* = "35.12": sta.lon* = "-85. 12": sta.alt* = "500.0" 

, Sample input (NAVSPASUR: Sat. 11111 data from SH0WALL, private comm. 

* Mr. Fred Lipp, July 1988): 
sat. id* = "11111" 

ecc* = "0.0005545": aincl* = "65.06057" 

peri* = "295.41470": anode* = "272.43497": manomaly* = "258.10682" 

motion* = "15.44194": decay* = "0" 

epochdate* = "1983.0201": epochtime* = "00.00" 

startdate* = "1983.0201": starttime* = "00.00" 

enddate* = "1983.0202": endtime* = "00.00" 

deltat* = "10" 

sta.id* = "001": sta.name* = "Daisy, Tenn." 

sta.lat* = "35.12": sta.lon* = "-85. 12": sta.alt* = "500.0" 

start : 

CLS 

PRINT "Select data input mode 1 , 2 , 3 or 4:" 

PRINT 

PRINT SPC(5); "Keyboard input for satellite k ground station:" 

PRINT SPC(IO); "1. Epoch of date." 

PRINT SPC(10); "2. NAVSPASUR Reference Date." 

PRINT 

PRINT SPC(5); "Data file input for satellite k ground station:" 

PRINT SPC(IO); "3. Epoch of date." 

PRINT SPC(10); "4. NAVSPASUR One-Line Charlie." 
data* = "" 

WHILE data* < "1" OR data* > "4" 
data* = INPUT* ( 1 ) 

WEND 

CLS 

PRINT "Select run mode 1 or 2:" 

PRINT 

PRINT SPC(5); "1. Create an ephemeris." 

PRINT SPC(5); "2. Find times satellite is above the horizon." 
case* = "" 

WHILE case* <> "1" AND case* <> "2" 
case* = INPUT* ( 1 ) 

WEND 

G0SUB input . s im . times 

SELECT CASE data* 

CASE "1", "2" 

CLS 

PRINT "Satellite Data:" 

PRINT 
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PRINT "Epoch date (yyyy.mradd) ? " ; epochdate* 

LOCATE 3, 29: INPUT v* : IF v* = "" THEN v* = epochdate* 
epochdate* = v*: CALL yrmoday(v*, yr k t mofc, dyfc) 

PRINT "Epoch time (hh.mmss) ? " ; epochtime* 

LOCATE 4, 29: INPUT v* : IF v* = "" THEN v$ = epochtime* 
epochtime* = v*: epochtime = decimalf(v$) 

CALL julian .day . number (mo* , dy*, yrJk, epochtime, jdepochO, epoch, tjepoch) 



LOCATE 6, 1: PRINT "Input mean elements of epoch:" 

PRINT "Satellite ID Number ? "; sat. id* 

LOCATE 7, 29: INPUT v* : IF v* = "" THEN v* = sat. id* 
sat. id* = v* 

PRINT "Eccentricity ? ecc* 

LOCATE 8, 29: INPUT v* : IF v* = "" THEN v* = ecc* 
ecc* = v*: ecc = VAL(v*) 

PRINT "Long, ascending node (deg) ? anode* 

LOCATE 9, 29: INPUT v*: IF v* = "" THEN v* = anode* 
anode* = v*: anode = VAL(v*) * rad 



PRINT "Inclination (deg) 7 "; aincl* 

LOCATE 10, 29: INPUT v$: IF v* = "" THEN v* = aincl* 
aincl* = v*: aincl = VAL(v*) * rad 



PRINT "Arg. of perigee (deg) 7 "; peri* 

LOCATE 11, 29: INPUT v*: IF v* = "" THEN v* = peri* 
peri* = v*: peri = VAL(v*) * rad 

PRINT "Mean anomaly (deg) ? "; manomaly* 

LOCATE 12, 29: INPUT v*: IF v* = "" THEN v* = manomaly* 
manomaly* = v*: manomaly = VAL(v*) * rad 



PRINT "Mean motion (rev/day) ? "; motion* 

LOCATE 13, 29: INPUT v*: IF v$ = "" THEN v* = motion* 

motion* = v* : motion = VAL(v$) * twopi / 1440# * convert to radians/minute 



PRINT "Decay (rev/day~2) 

LOCATE 14, 29: INPUT v*: IF v* = 

1 convert to radians/minute * 2 : 
decay* = v*: decay = VAL(v*) * tuopi / (1440#) * 2 



"; decay* 

THEN v* = decay* 



LOCATE 16, 1: PRINT "Press any Xey to continue." 
FOR i 7. = 1 TO 10: c* = INKEY*: NEXT 
c* = INPUT$( l) 



CLS 

LOCATE 1, 1: PRINT "Station Coordinates: 
PRINT 



PRINT "Station Identification Number 
LOCATE 3, 42: INPUT v*: IF v* = "" THEN v* 
sta.id* = v* 



? " ; sta . id* 
sta . id* 



PRINT "Station Location or Name 

LOCATE 4, 42: INPUT v* : IF v* = "" THEN v* 

sta. name* = v*: sta.lat = VAL(v*) * rad 



? "; sta. name* 
sta . name* 



PRINT "Station latitude (dd.mmss, South minus) ? sta.lat* 
LOCATE 5, 42: INPUT v* : IF v* = "" THEN v* = sta.lat* 
sta.lat* = v*: sta.lat = VAL(v*) * rad 



PRINT "Station longitude (ddd.mmss, Vest minus) ? "; sta. Ion* 
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LOCATE 6, 42: INPUT v$ : IF v$ = "" THEN v$ = sta.lon$ 
sta.lon$ = v$: sta.lon = VAL(v$) * rad 

PRINT "Station altitude (feet) ? sta.alt$ 

LOCATE 7, 42: INPUT v$ : IF v$ = "" THEN v$ = sta.alt* 
sta.alt* = v$: sta.alt = VAL(v*) 

LOCATE 9, 1: PRINT "Press any key to continue." 

FOR it = 1 TO 10: c* = INKEY* : NEXT 
c* = INPUT$( 1 ) 

OPEN "output . dat " FOR APPEND AS #3 
GOSUB compute 
CASE "3" 

OPEN "station.dat" FOR INPUT AS #1 
OPEN "output . dat" FOR APPEND AS #3 

DO UNTIL EOF(l) 

INPUT #1, sta. id*, sta.lat*, sta.lon*, sta.alt*, sta.narae$ 
sta.lat = VAL ( sta . lat*) * rad 
sta.lon = VAL(sta . lon$) * rad 
sta.alt = VAL(sta . alt$) 

OPEN "elements . dat" FOR INPUT AS #2 

DO UNTIL E0F(2) 

INPUT #2, sat. id*, epochdate*, epochtime*, ecc*, anode*, aincl*, _ 
peri*, manomaly*, motion*, decay$ 

CALL yrmoday(epochdate* , yr& , mo£ , dy&) 
epochtime = decimal# ( epochtime*) 

CALL julian . day . nuraber(mo£ , dy b t yr&, epochtime, jdepoch(), epoch, 
t j epoch) 
ecc = VAL(ecc$) 
anode = VAL(anode*) * rad 
aincl = VAL(aincl$) * rad 
peri = VAL(peri$) * rad 
manomaly = VAL(manomaly$) * rad 

motion = VAL(motion*) * twopi / 1440#' convert to radians/minute 
decay = VAL(decay$) * twopi / (1440#) “ 2 1 convert to rad/mm“2 

jd(l) = jdstart(l) 
jd(2) = jdstart(2) 
jd = jd( 1) ♦ jd(2) 
time. of. day = starttime 

GOSUB compute 

LOOP 
CLOSE 2 
LOOP 

CASE "4" 

OPEN "station.dat" FOR INPUT AS #1 
OPEN "output .dat" FOR APPEND AS #3 

DO UNTIL EOF(l) 

INPUT #1, sta . id$ , sta . lat$ , sta.lon*, sta.alt*, sta.name$ 
sta.lat = VAL(sta . lat$) * rad 
sta.lon = VAL(sta . Ion*) * rad 
sta.alt = VAL(sta . alt*) 

OPEN "elements. nss" FOR INPUT AS #2 
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DO UNTIL EOF( 2) 

INPUT #2, c$ 

sat.idS = MID$(c$, 1, 5) 
manomalyS = MIDS(cS, 6, 9) 
motionS = MIDS(cS, 14, 9) 
decayS = MID$(cS, 22, 7) 
eccS = MIDS(c$, 28, 9) 
periS = MIDSCcS, 36, 9) 
anodeS = MID$(cS, 44, 9) 
ainclS = MIDSCcS, 52, 9) 

epochdateS = "19" ♦ MIDSCcS, 60, 2) ♦ ♦ MIDSCcS, 62, 4) 

epochtimeS = "00.00" 

CALL yrmoday C epochdateS , yr it, mot, dyt) 
epochtime = decimal# CopochtimeS) 

CALL Julian . day . number (mo* , dyt, yrt, epochtime, jdepoch(), epoch, 
t j epoch) 
ecc = VAL(eccS) 
anode = VALCanodeS) * twopi 
aincl = VALCainclS) * twopi 
peri = VAL(periS) * twopi 
manomaly = VALCmanomalyS) * twopi 

motion = VALCmotionS) * Herg.min ' convert to radians/minut e 

decay = VALCdecayS) 

IF decay > .5# THEN decay = decay - 1# 

decay = decay * .00001 * Herg.min * 2 ' convert to rad/min~2 

j d C 1 ) = jdstart(l) 
j d C 2 ) = jdstart(2) 
jd = jd C 1 ) + j d C 2 ) 
time. of. day = starttime 

G0SUB compute 

LOOP 
CLOSE 2 
LOOP 

END SELECT 
CLOSE 

END 



> * 



* 



* 



♦ 



compute : 

1 initialize times, 
tj = tjstart 

begin. time = Cstartdate - epoch) * 1440# 
end. time = Conddate - epoch) * 1440# 



'julian centuries from J2000.0 
'minutes 
' minutes 



elem( 1 ) 


= 


amaj or 


elem(2) 


= 


ecc 


elem(3) 


= 


tzero 


elem(4) 


= 


aincl 


elem(5) 


= 


anode 


elem(6) 


= 


peri 


elem(7) 


= 


motion 


elem(8) 


= 


manomaly 


elem( 9) 


= 


decay 


SELECT CASE < 


dat 


;aS 



'semi-major axis 
'eccentricity 

'time of perifocal passage 
'inclination (radians) 

'longitude of ascending node (radians) 
'argument of perifocus (radians) 

'mean motion (revolut ions/day ) 

'mean anomaly (radians) 

'decay constant (radians/minute" 2) 



CASE "1" , "3" 

' Since the earth is treated as non-rotating, elem(5) must be corrected for 
' the accrued rotation from epoch to J2000.0 
elem(5) = elem(5) - 15# * gmst(tj epoch , 0#) * rad 
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CASE "2", "4" 

’ NAVSPASUR one-line Charlie correction 

CALL julian .day .number(l , 1, 1985, 0, jd . Charlie () , dummy, t j . Charlie) 
vetc = (tjepoch - tj. charlie) * 36525# - (l# - gms t(tj . charlie , 0#) / 24#) 
* .9972695664# 

vetc = dmod(wetc * 1440# * we, twopi) 
elem(5) = elera(5) + vetc 

elera(5) = elem(5) - 15# * gmst(tj epoch , 0#) * rad 
END SELECT 



CALL sta.p.coord(8ta.lat , sta.lon, sta.alt, sta.xQ) 

FOR i'/. = 1 TO 3 

sta.z(iX) = sta.x(iy, , 2) 
sta.R(iX) = sta.x(i%, 1) 

NEXT 

CALL pos .updat e( 1 , elem(), begin. time, sat.xQ, sat.xdQ, elemuQ, r, ea) 
CALL pos . update(2 , elemQ, begin. time, sat.xQ, sat.xd(), elemu(), r, ea) 

CLS 

PRINT "Satellite: sat.id$; " Station: sta.id$; " sta.nameS 

PRINT 

SELECT CASE case$ 



CASE "1" 

tm = starttime ’ UT in hours 

CALL apo .peri .gee(elem() , elemuO) 

CALL jd . to .date( jd( ) , mk , d& , y&, h# , mo$, day$) 

PRINT 

PRINT day$; dk; mo$ ; yk 

CALL output4(0, yr&, mo& , da&, time. of. day, sat.lt, sat. In, sat.ht, el, az , 
rng , langle , head, sta.id$, sat.idS) 



* time-loop 

FOR time = begin. time TO end. time STEP delta. t 

CALL pos . update (2 , elem(), time, sat.xQ, sat.xd(), elerau(), r, ea) 

’compute spherical coordinates 

CALL sat . sphere . coord(sat . x () , r, t j , tm, sat.dc, sat. In, sat.ra) 

’satellite heading 

CALL sat.head(r, sat.x(), sat.xdQ, head) 

’earth ground track 

CALL ground .track(r, sat.dc, sat.lt, sat.ht) 

’compute station-satellite relations 

CALL sta . sat(sta.x() , sat.xQ, rng, az , el, langle) 

rng = rng * ae 

CALL output4(l, yk , m&, d&, time. of. day, sat.lt, sat. In, sat.ht, el, az , rng, _ 
langle, head, sta.id$, sat.id$) 

’update time 

CALL t ime .update(0 , delta. t. hr, tm, jd, jd(), t j , time. of. day, m& , dk , y &, _ 
h#, mo$, day$) 



NEXT time 
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PRINT : PRINT "Press any key to continue." 
FOR iX = 1 TO 10: c$ = INKEYS: NEXT 
c$ = INPUTS ( 1 ) 



CASE "2" 



CALL output4(0, yt , m& , dfc, time. of. day, eat.lt, sat. In, sat.ht, el, az , rng , _ 
langle, head, sta.id$, sat.idS) 

tm = 0 

period = tvopi / motion 
ea = 0 



time = begin. time 
initialize : 

CALL pos . update(2 , elemQ , time, sat.x(), sat.xd(), elemu(), r, ea) 

CALL culmnate (sta.z() , we, time, ea, coseta, elemuQ) 

IF coseta < 0 THEN 
ea = ea + pi 

CALL culmnate(sta . z( ) , we, time, ea, coseta, elemuO) 

END IF 

time = (ea - elem(2) * SIN(ea)) / motion + elemu(3) 
dt.hr = (time - begin. time) / sixty 

CALL time . update(0 , dt.hr, tm, jd, jd(), tj , time. of. day, m& , dt , y t, h# , _ 
mo$, day$) 

IF time > begin. time ♦ period THEN 
time = time - period 

CALL time .updated , -period, tm, jd, jd(), t j , time. of. day, m& , dt , yt , _ 
hf , mo$ , day$) 

GOTO initialize 
END IF 

again : 

FOR iX = 1 TO 2 

CALL pos . update (2 , elem() , time, sat.x(), sat.xdO, elemuO, r, ea) 

CALL culmnate(8ta .z() , we, time, ea, coseta, elemuO) 
dt.min = (ea - elem(2) * SIN(ea)) / motion ♦ elemu(3) - time 
time = time ♦ dt.min 

CALL time . update( 1 , dt.min, tm, jd, jd(), t j , time. of. day, mt , d& , y&, _ 
hf, mo$ , day$) 

NEXT 

time.cul = time 



CALL sta . sat ( sta . x ( ) , sat.x(), rng, az, el, langle) 

IF el > 0 AND time >= begin. time THEN 

CALL pos .update(2, elem(), time, sat.xQ, sat.xdO, elemuO, r, ea) 
ea = ea - .335 
FOR iX = 1 TO 2 

CALL ri se . 8 et (sta . z( ) , sta.R(), we, time, ea, elemuO) 
dt.min = (ea - elem(2) ♦ SIN(ea)) / motion ♦ elemu(3) - time 
time = time ♦ dt.min 

CALL time . update( 1 , dt.min, tm, jd, jd() , tj, time. of. day, mt , dt , _ 
yt , hf , mo$, day$) 

CALL pos .update(2 , elem() , time, sat.xO, sat.xdO, elemuO, r, ea) 

NEXT 

time. rise = time 

CALL 8 ta . sat ( sta . x () , sat.xO, rng, az , el, langle) 

CALL sat .sphere. coord(sat.x() , r, t j , tm, sat.dc, sat. In, sat.ra) 

CALL ground . track(r , sat.dc, sat.lt, sat.ht) 

CALL sat . head(r , sat.xO, sat.xdO, head) 

CALL jd . t o . dat e ( jd( ) , mt , dt , yt , hf , mo$ , day$) 

CALL output4 ( 1 , yt, mt , dt , time. of. day, sat.lt, sat. In, sat.ht, el, _ 
az , rng * ae, langle, head, sta.id$, sat.idS) 
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time.incr = dmod(time . cul - time. rise, delta. t) 
time = time ♦ time.incr 

CALL time . updat e(l , time.incr, tm, jd, jd(), t j , time, of. day, mb, db , _ 
yb t ht, mo$, day$) 

CALL pos . update(2 , elemQ, time, sat.xO, sat.xdO, elemu(), r, ea) 

CALL sta. sat(sta.x() , sat.x(), rng , a z, el, langle) 

CALL sat . sphere . coord(sat ,x() , r, t j , tm, sat.dc, sat. In, sat.ra) 

CALL ground . track(r , sat.dc, sat.lt, sat.ht) 

CALL sat.head(r, sat.xC), sat.xdO, head) 

CALL jd.to.date( jd() , mk, d b t yk, hi, mo$, day$) 

CALL output4(l, y b, mb , d& , time. of. day, sat.lt, sat. In, sat.ht, el, _ 
az, rng * ae , langle, head, sta.id$, sat.id$) 



DO WHILE el > 0 

time = time ♦ delta. t 

CALL time . update(0 , delta. t. hr, tm, jd, jd(), t j , time. of. day, mb , db , _ 
y&, ht, mo$, day$) 

CALL pos .update(2 , elem(), time, sat.x(), sat.xdO, elemu(), r, ea) 

CALL sta . sat (sta. x() , sat.xO, rng, az, el, langle) 

IF el < 0 THEN 
EXIT DO 
END IF 



LOOP 



CALL sat . sphere . coord(sat . x() , r, t j , tm, sat.dc, sat. In, sat.ra) 
CALL ground . track(r, sat.dc, sat.lt, sat.ht) 

CALL sat.head(r, sat.xO, sat.xdO, head) 

CALL jd . to . dat e( jd() , m& , db , yb , hf, mo$ , day$) 

CALL output4(l, yb, mb, d& , time. of. day, sat.lt, sat. In, sat.ht, el, 
az , rng * ae , langle, head, sta.id$, sat.id$) 



save. time = time 
FOR 1 7. = 1 TO 2 

CALL pos . update (2 , elem() , time, sat.xO, sat.xdO, elemuO, r, ea) 

CALL rise . set (sta. z() , sta.R(), we, time, ea, elemuO) 
dt.min = (ea - elera(2) * SIN(ea)) / motion ♦ elemu(3) - time 
time = time ♦ dt.min 

CALL time . updat e( 1 , dt.min, tm, jd, jd(), t j , time. of. day, mb, d&, yb, _ 
ht, mo$, day$) 

NEXT 

CALL pos . update (2 , elem(), time, sat.xO, sat.xdO, elemuO, r, ea) 

CALL sta. sat(sta.x() , sat.xO, rng, az , el, langle) 

CALL sat . sphere . coord(sat . x() , r, t j , tm, sat.dc, sat. In, sat.ra) 

CALL ground . track(r, sat.dc, sat.lt, sat.ht) 

CALL sat.head(r, sat.xO, sat.xdO, head) 

CALL jd . to .date( jd() , mb , d&, yb, ht, mo$ , day$) 

CALL output4(l, yb, m& , db , time. of. day, sat.lt, sat. In, sat.ht, el, _ 
az , rng * ae, langle, head, sta.id$, sat.id$) 

PRINT 
END IF 

time = time. cul ♦ period 

dt.hr = (time * begin. time) / sixty - tm 

CALL time .update(0 , dt.hr, tm, jd, jd(), t j , time. of. day, m& , d&, y&, ht , . 
mo$ , day$) 

IF time < end. time THEN GOTO again 



END SELECT 
RETURN 

input . sim. times : 

CLS 

PRINT "Simulation-Time Data:" 

PRINT 

PRINT "Starting date (yyyy .mmdd) ? startdate$ 
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LOCATE 3, 29: INPUT v$ : IF v$ = THEN v$ = etartdatel 
startdate$ = v$: CALL yrmoday(v$, yr *, mo*, dy*) 

PRINT "Starting time (hh.mmss) ? starttimeS 

LOCATE 4, 29: INPUT v$ : IF v$ = "" THEN v$ = starttimeS 
starttimeS = v$ : starttime = decimalf(v$) 

CALL julian. day .number (molt, dy*, yr*, starttime, jdstartO, startdate, 

tj start ) 

j d ( 1 ) = jdstart(l) 
jd(2) = jdstart (2) 
jd = jd(l) ♦ jd(2) 
time. of. day = starttime 

PRINT "Ending date (yyyy.mmdd) ? enddateS 

LOCATE 5, 29: INPUT vS : IF vS = " M THEN vS = enddateS 
enddateS = vS: CALL yrmoday(v$, yr* , mo*, dy*) 

PRINT "Ending time (hh.mmss) ? endtimeS 

LOCATE 6, 29: INPUT vS : IF vS = "" THEN vS = endtimeS 
endtimeS = vS : endtime = decimalt(vS) 

CALL julian .day .number (mo* , dy&, yr*, endtime, jdend(), enddate, tjend) 

PRINT "Time increment (minutes) ? deltatS 

LOCATE 7, 29: INPUT vS : IF vS = "" THEN vS = deltatS 

deltatS = vS : delta. t = VAL(v$): delta. t. hr = delta. t / sixty 

LOCATE 9, 1: PRINT "Press any key to continue." 

FOR i*/. = 1 TO 10: c$ = INKEYS: NEXT 
cS = INPUT$(1) 



RETURN 



K*»*«x»xa 



>*********************************************X******X*****X* 

FUNCTION acos# (xf) ' 02-10-88 Rev. 02-24-88 

’ The arccosine function derived from the ATN function with no singularities. 
' The angle is returned in the (0,pi) interval. 

'Note that in relational comparisons, -1 is "true" and 0 is "false". 

CONST pi = 3.1415926535897931, eps = ID-33 
IF ABS ( x#) <= 1# THEN 

acosl = ATN(SQR(l# - xf * xf) / (x# - eps * (xf = Of))) - pi * (x# < Of) 
ELSE 

PRINT 

PRINT "Error in acos function. ABS(arg) > 1" 

PRINT "arg = " ; xf 
PRINT 
END IF 

END FUNCTION 



>******♦***********★***♦********★******************************************** 
SUB apo. peri. gee (elemi(), elemufO) ' 09-14-88. Rev 09-14-88. 

' Compute location and height of apogee and perigee above oblate earth. 

1 NOTE: (1) The output from this routine is approximately valid for the 
' current orbit only. In particular, the longitude is not corrected 

* for the earth's rotation during the current orbit. Also, the 

' ascending node and argument of perigee are not corrected for 

' their rates of change during the current orbit. 

' (2) For nearly circular orbits, the height of the satellite will 

' probably get lower than the perigee value and higher than the 

' apogee value at locations other than the perigee and apogee 

' positions this is an effect of the earth's oblateness. 
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CONST pi = 3.141592653589793#, rad = pi / 180#, twopi = pi ♦ pi 



* spherical earth values: 

apogee# = elemu#(l) * (1# - elem#(2)) 

perigee# = elemu#(l) * (1# ♦ elem#(2)) 

' longitudes and declinations: 
sinci# = SIN( elem#(4)) 
speri# = SIN( elemu#(6)) 
apo.dc# s asin#(speri# * sinci#) 

apo.ln# = elemu#(5) ♦ atan2#(C0S(elem#(4) ) * speri#, C0S(elemu#(6) ) ) 
apo.ln# = dmod#(apo .In#, twopi) 

IF apo.ln# > pi THEN apo.ln# = apo.ln# - twopi 

' next are approximations should update by peri- k nodal-rates 
* for 1/2 revolution: 

peri.dc# = -apo.dc# 

peri. In# = dmod(apo . In# + pi, twopi) 

IF peri. In# > pi THEN peri. In# = peri. In# - twopi 

CALL ground. track(apogee# , apo.dc#, apo.lt#, apo.ht#) 

CALL ground . track(perigee# , peri.dc#, peri. It#, peri.ht#) 



PRINT 

PRINT "Apogee : " 

PRINT USING "k ####.# *" ; 
PRINT USING "k ###.####"; 
PRINT USING "k ###.####"; 
PRINT 

PRINT "Perigee:" 

PRINT USING "k ####.# k " ; 
PRINT USING 11 k ###.####"; 
PRINT USING "* ###.####"; 
PRINT 



"Height = 
"Latitude = 
"Longitude = 



"Height = 
"Latitude = 
"Longitude = 



apo.ht#; " km. " 

" ; apo.lt# / rad 
" ; dmod(apo.ln# / rad, 360#) 



peri.ht#; " km. " 

"; peri. It# / rad 
"; dmod(peri . In# / rad, 360#) 



IF apo.ht# < 0 THEN 

PRINT "WARNING: Satellite is suborbital and will impact earth." 
PRINT "Do you want to continue? Yes or No." 

PRINT 

WHILE c$ <> "y" AND c$ <> "n" 
c$ = LCASE$ ( INPUTS ( 1 ) ) 

WEND 

IF c$ = "n" THEN END 
END IF 
END SUB 



FUNCTION asin# (xt) * 02-10-88 Rev. 02-24-88 

* The arcsine function derived from the ATN function with no singularities. 

* The angle is returned in the (-pi, pi) interval. 

'Note that in relational comparisons, -1 is "true" and 0 is "false". 

CONST eps = ID-33 
IF ABS(x#) <= 1# THEN 

asin# = ATN(x# / (SQR(l# - x# * x#) - eps * (ABS(x#) = 1#))) 

ELSE 

PRINT 

PRINT "Error in asin function. ABS(arg) > 1" 

PRINT "arg = " ; x# 

PRINT 
END IF 

END FUNCTION 

FUNCTION atan2i (yf, x#) ’ 02-10-88 Rev. 02-24-88 
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' The two argument quadrant determining arctangent function. 

* The angle is returned in the (-pi, pi; interval. 

* Note that in relational comparisons, -1 is "true” and 0 is “false". 
CONST pi = 3.141592653589793#, eps = ID-33 

atan2# = ATN(y# / (x# - eps * (x# = 0#))) - pi * (x# < 0#) * (SGN(yf) _ 

- (y# = 0#)) 



END FUNCTION 



’***************************************+*+*****************************+****** 

SUB culmnate (zv#(), wet, t0#, eat, cosetal, elemu#()) 



Compute eccentric anomaly for the time of culmination of a satellite 
01-30-89. Rev. 02-02-89 C 1000. 



NOTE: This particular form of the culminate function was written from the 

development of 30 Jan 89. It assumes a non-rotating earth so that the 
station longitude is fixed. Rotation is accomplished by subtracting the 
accumulated rotation from the right ascension of the ascending node of 
the satellite. This is the method used to compute the ephemens and is 
also equivalent to the method used by NAVSPASUR. 

lambda - station geodetic longitude 

phi = station geodetic latitude 

we = earth's rotation (sidereal rate of change) 

tO = time the elements, elemuQ, were last updated 

capt = T -- time of latest perifocal passage 

mot - mean motion 

ea = eccentric anomaly E 

ecc = eccentricity e 

amjr = semimajor axis (cancels out in all formulas, so not included) 



CONST pi = 3.141592653589793#, rad = pi / 180#, twopi = pi ♦ pi 



ecc# = elemu#(2) 
capt# = elemu#(3) 
incl# = elemu#(4) 
node# = elemu#(5) 
peri# = elemu#(6) 
mot# - elemu#(7) 



DIM pv# ( 3) , pvd#(3) , pvdd#(3), qv#(3) , qvd#(3) , qvdd#(3) 

DIM rvec#(3), rvecd#(3), rvecdd#(3) 

' zv = Z-vector, unit normal at the station 

' pv = P-vector, qv = Q-vector (first two columns of the Euler matrix) 

' rvec = R-vector (satellite position), rvecd = dR/dE, rvecdd = d2R/dE2 

' rv = r magnitude of the R-vector, rvd = dr/dE, rvdd = d2r/dE2 

* gamma = node - accrued rotational motion 

cosw# = C0S(peri#) 
sinw# = SIN(peri#) 
cosi# = C0S(inc 1#) 
sini# = SIN(incl#) 



corr# = 1 

WHILE ABS( corr#) > .0001 
cosea# = C0S(ea#) 
ecosea# - ecc# * cosea# 
sinea# = SIN(ea#) 
esinea# = ecc# * sinea# 

gamma# = node# - we# * ((ea# - esinea#) / mot# ♦ capt# - t0#) 
gammad# = -we# * (1# - ecosea#) / mot# 
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gammadd# = -we# * esineaf / mot# 

cosg# = C0S(gamma#) 
sing# = SIN(gammat) 

rtecc# = SQR(1# * ecc# * ecct) 

rv# = 1# - ecosea# 
rvd# = esinea# 
rvdd# = ©cosea# 



xwt = coseaf - ecc# 
xwd# = -sinea# 
xwdd# = -cosea# 
yw# = rtecc# * sin©a# 
ywd# = rtecc# * cosea# 
ywdd# = -y u# 

pv#(l) = cosw# * cosg# - sinw# * sing# * cosi# 
pv#(2) = cos u# * sing# + sinw# * co sg# * cos i# 
pv#(3) = sinw# * sini# 

qv#( 1) = -sinw# * cosg# - cosw# * sing# * cosi# 
qv#(2) = -sinw# * sing# + cosw# * cosg# * cosi# 
qv#(3) = cosw# * sini# 

pvd#(l) = -pv#( 2) * gammad# 
pvd#(2) = pv#(l) * gammad# 
pvd# (3) = 0# 

qvd# ( 1 ) = -qv#(2) * gammad# 
qvd#(2) = qv#(l) * gammad# 
pvd#(3) = 0# 

pvdd#(l) = -pvd#(2) * gammad# - pv#(2) * gammadd# 
pvdd#(2) = pvd#(l) * gammad# + pv#(2) * gammadd# 
pvdd#(3) = 0# 

qvdd#(l) = -qvd#(2) * gammad# - qv#(2) * gammadd# 
qvdd#(2) = qvd#(l) * gammad# + qv#(2) * gammadd# 
pvdd#(3) = 0# 



FOR iX = 1 TO 3 



rvec#(i'/.) = xw# * pv#(i‘/*) + yw# * qv#(iX) 

rv©cd#(iX) = xwd# * pv#(iX) * xw# * pvd#(iX) + ywd# * qv#(i'/.) _ 

♦ yw# * qvd#(iX) 

temp# = xwdd# * pv#(iX) ♦ 2# * xwd# * pvd#(i50 + xw# * pvdd#(i%) 
rvecdd#(i'/.) = t©mp# + ywdd# * qv#(i>0 + 2# * ywd# * qvd#(i'/.) _ 

♦ yw# * qvdd#(iX) 



NEXT 



cos©ta# = dot#(zv#(), rv©c#()) / rv# 

fe# = dot#(zv#(), rvecd#Q) - rvd# * coseta# 

rrd# = rvd# / rv# 

fed# = dot#(zv#(), rvecdd#()) + coseta# * (rrd# - rvdd#) 
fed# = fed# - rrd# * dot#(zv#(), rvecd#()) 

corr# = fe# / fed# 

IF corr# > pi THEN 

ea# = ea# ♦ twopi 

ELSE 

ea# = ea# - corr# 

IF coseta# < 0 THEN ea# = ea# ♦ pi 
END IF 



WEND 
END SUB 
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> **************************************************************************** 
FUNCTION decimal# (v$) 1 03-21-88. Rev 03-21-88. 

’ Convert ddd.mmssff or hh.mmBsfff to decimal degrees or decimal hours. 

ii = INSTR(v$, 

IF ix = 0 THEN 

decimal# = VAL(v$) 

ELSE 

x# = VAL(LEFT$(v$, ix)) 
sn = 0 

IF x# < 0# THEN 
sn = -1 
x# = -x# 

END IF 

w$ = v$ ♦ "0000" 

y# = VAL(MID$(w$, ix ♦ 1, 2)) 

z# = VAL(MID$(w$, ix ♦ 3, 2) ♦ " . M ♦ RIGHT$(w$, LEN(w$) - ix - 4)) 
x# = (z# / 60# ♦ y#) / 60# ♦ x# 

IF sn THEN x# = -x# 
decimal# = x# 

END IF 

END FUNCTION 

***************************************************************************** 
FUNCTION dmod# (x#, m#) » 02-10-88 Rev. 02-24-88 

’ Provides x MOD m for real-valued functions. 

IF m# <> 0# THEN 

dmod# = x# - m# * INT(x# / m#) 

ELSE 

PRINT 

PRINT "Error in dmod function. Modulus cannot be zero." 

PRINT 
END IF 

END FUNCTION 

***************************************************************************** 
FUNCTION dot# (a#(), b#()) 

’ Dot product of two vectors . 

dot# = a#(l) * b#( 1 ) ♦ a# (2) » b*(2) + a#(3) • b# ( 3 ) 

END FUNCTION 

'*********************************************************************»****** 
SUB euler (peri#, node#, incl#, psn#(), vel#()) 

* 09-05-88. Rev. 01-16-89. 

’Compute the first two columns (p k q) of the Euler matrix. 

’Rotate in-plane position and velocity components to rectangular equatorial 
’ coordinates. 

’psn() and vel() are both input and output vectors. It is assumed that the 
’ third component of both psn() and vel() is zero on input. 



cw# 


= 


C0S(peri#) 






BW# 


= 


SIN(peri#) 






cn# 


= 


C0S(node#) 






sn# 


= 


SIN(node#) 






ci# 


= 


C0S( incl#) 






si# 


= 


SIN(inclf) 






px# 


= 


cw# * cn# - 


sw# * sn# * 


ci# 


py# 


= 


cw# * sn# ♦ 


sw# * cn# * 


Cl# 


pz# 


= 


sw# * si# 






qx# 


= 


-BW# * cn# 


- cw# * sn# * 


Cl# 


qy# 


= 


-sw# * sn# 


♦ cw# * cn# * 


ci# 


qz> 


= 


cw# * si# 
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tit = psn#( 1 ) 
t2t = psn#(2) 

psnt(l) = pxt * tit + qxt * t2t 

psnt(2) = pyt * tit ♦ qyt * t2t 

psnt(3) = pzt * tit + qzt * t2t 

tit = velt(l) 
t2t = velt (2) 

velt(l) = pxt * tit + qxt * t2t 

velt(2) = pyt * tit + qyt * t2t 

velt(3) = pzt * tit + qzt * t2t 

END SUB 

t ****************************************************************************** 

FUNCTION gmstt (tut, tmt) 1 03-03-88. Rev 03-03-88. 

' Compute Greenwich mean sidereal time. 

' tut = number of centuries of 36525 days of universal time elapsed since 
' 2000 Jan 1, 12h UT1 (JD 2451545 UT1). 

' tmt = time of day 

' The Astronomical Almanac, 1984, pp S13-S15. 

1 Almanac for Computers 1988, pp B2-B3. 

CONST pi = 3. 14 15926535897 93 t, rad = pi / 180t 

CONST cO = 24110.54841#, cl = 8640184 . 812866t , c2 = 9 . 310400000000001D-02 
CONST c3 = -.0000062#, hourpersec = It / 3600t, meantoapp = -.00029# 

CONST dO = 125.004452#, dl = -.0529538#, d2 = .002071# 

tt = C C ( c 3 * tut + c2) * tut + cl) * tut + cO) * hourpersec + tmt' GMST 
'The following two lines will convert GMST to CAST 
'lunarnodet = (d2 * tut + dl) * tut + tOt 
'tt = tt + meantoapp * SIN (lunarnodet * rad) 
gmstt = dmodt(tt, 24t) 

END FUNCTION 

******************************************************************************* 
SUB ground. track (rt, sat.dct, sat.lt t, sat.htt) '09-14-88. Rev 09-14-88. 

CONST ae = 6378.14# ' earth's equatorial radius (km.) 

CONST fl = It / 298.257# ' earth's flattening factor 

CONST e.sqr = (2t - fl) * fl ' earth's ellipticity squared 

CONST f 2 = (It - fl) * (It - fl) 

rkmt = rt * ae 
pst = sat .dct 
el: cpt = COS(pst) 

ret = ae * SQR((lt - e.sqr) / (it - e.sqr * cpt * cpt)) 
rert = ret / rkmt 

sat.lt t = ATN(TAN(pst) / f2) ' satellite latitude 

snt = SIN(sat.lt# - pst) 

hr# = SQR(lt - rert * rert * snt * snt) - rert * COS(sat.ltt - pst) 
pnt = sat.dct - asin(hrt * snt) 

IF ( ABS(pnt - pst) > .0000001) THEN pst = pnt: GOTO el 
sat.htt = hr# * rkmt ' satellite height 

END SUB 



******************************************************************************* 



SUB hr. min. sec (xt, hrX, min?., 
' Convert decimal hours to: 



sect, nX) 
hh mm ss 
hh mm ss.f 
hh mm ss . f f 
hh mm ss.fff 



if nX = 0 , 
if nX = 1, 
if nX = 2, 
if nX = 3. 



CONST tO = It, tl = 10#, t2 = 100#, t3 = 1000# 
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CONST rO = If / 7200#, rl = rO / tl, r2 = r 1 / t2 , r3 = r2 / t3 

SELECT CASE nX 
CASE 0 

round - rO 
trunc = tO 
CASE 1 

round = rl 
trunc = tl 
CASE 2 

round = r2 
trunc = t2 
CASE 3 

round = r3 
trunc = t3 
CASE ELSE 
PRINT 

PRINT "Error in hr, min. sec subroutine. Must have 0 <= nX <= 3." 
PRINT "nX = nX 
PRINT 
EXIT SUB 

END SELECT 

vf = ABS(xf) ♦ round 
hrX = INT(vf) 
vf = (vf - hrX) * 60# 
minX = INT(vf) 

sec# = INT((60# * (v# - minX) * trunc)) / trunc 
END SUB 

J*********************************************************************xx*»** 

SUB jd. to. date (jd#(), mft , dt , yt. h# , mo$ , day$) > 02-24-88 Rev. 03-02-88 
’ Convert Julian Day Number to month number, day, year, hour, month name and 
’ day of the week. Limited to the Gregorian Calendar. 

’ Fliegel it VanFlandern, Comm. ACM, Vol. 11, No. 10, Oct. 1968, pg . 657. 

CONST «$ = "MonTueWedThuFriSatSun" 

CONST dd$ = " JanFebMarAprMayJunJulAugSepOctNovDec" 

j 4 Ac = INT(jd#(l) + jd#(2) + .5) 

h# = (( jd#(l) - j 4 Ac ) + jd#(2) + .5) « 24# 

lft = INT ( j4t + 68569) 

nft = INT(4 * lft / 146097) 

lft = lft - INT ( (146097 * nft + 3) / 4) 

yft = INT (4000 * (lft + 1) / 1461001) 

It = It - INT ( 1461 * yft / 4) ♦ 31 

mft = INT(80 * lft / 2447) 

dft = It - INT(2447 * mft / 80) 

It = INT (mft / 11) 

mft = mft + 2 - 12 * It 

yft = 100 * (nft - 49) + yft + lft 

vny. = j4ft - 7 « INT(j4t / 7) + 1 
day$ = MID$( w$ , 3 * un'/. - 2, 3) 
mo$ = MID$(dd$ , 3 * mft - 2, 3) 

END SUB 

l*M**M**M*******************#****M****M*****M************************* 

SUB julian .day . number (m&, dk , y & , utf, jdf(), djf, tjf) 

» 02-25-88 Rev. 03-23-88 

} Conversion of calendar date: m* = month, df = day, yk = 4-digit year, 

’ and utf = Universal Time, to Julian Day Number. Year must be between 1801 
1 and 2099. 
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’ Almanac for Computers, 1988, pg. B-2. 

’ Note that jd# is a subscripted double precision array. At any time, the 
* Julian Day Number is given by jd#(l) ♦ jd#(2). 

’ For ease of programming, the user may put the entire epoch in jdf(l) and set 
1 jdf (2) = Of. 

’ For maximum accuracy, set jdf(l) = the moat recent midnight at or before the 
’ epoch and set jd#(2) = the fractional part of a day elapsed between jdf(l) 

’ and the epoch. As an example, compute the number of days from epoch 

’ J2000.0 = JD 2451545.0 using the code: 

1 days# = ( jd#( 1) - 2451545.0#) ♦ jd#(2) 

’ Note, particularly, the extra pair of parenthesis. 

’ dj# = days and fraction from J2000.0. tj# = julian centuries from J2000.0 



CONST cl = 1721013.5#, c2 = 190002.5#, jcent = 36525# 

IF (1801 <= yb) AND (y k <= 2099) THEN 

jd#( 1 ) = 367 * yk - INT(7 * (y k * INT((mfc ♦ 9) / 12)) / 4) . 

♦ INT( (275 * mb) / 9) ♦ db ♦ cl - .5# * SGN(100 * yb ♦ mb - c2) ♦ .5# 
jd#(2) = ut# / 24 
dj# = ( jd#( l) - 2451545#) + jd#(2) 
tj# = dj# / jcent 
ELSE 

PRINT 

PRINT "Error in j ulian . day . number . Year is not between 1801 and 2099." 

PRINT "Year = "; y b 
PRINT 
END IF 
END SUB 

tit********************************************************************** ******+ 

SUB kepler (m#, eel, ea#) ' 02-25-88. Rev 02-26-88 

’High precision, non-iterative, solution to Kepler’s equation for the ellipse. 
’Accuracy is 10E-18 using all terms or 10E-15 if the last term is omitted. 
’Although not as compact as the Newton-Raphson procedure, one square root, one 
’cube root, and only two trigonometric functions are evaluated thus malting this 
’procedure extremely fast. Can be extended to the hyperbolic case - see Ref. 
’Ref.: S. Mikkalo , Cel. Mech., Vol. 40, 1987, pp 329-334. 

’m# = mean anomaly, ec# = eccentricity and ea# = eccentric anomaly. 

CONST pi = 3.141592653589793#, twopi = pi ♦ pi , third = 1# / 3# 

CONST r3 = 1# / 6#, r4 = 1# / 24#, r5 = 1# / 120# 

IF ec# >= 1# OR ec# < 0# THEN 
PRINT 

PRINT "Error in Kepler. The eccentricity is not in the [0,1) interval." 
PRINT "Eccentricity = ec# 

PRINT 

ELSE 

m# = m# - twopi * INT(m# / twopi)’ assure that m# is 

IF m# > pi THEN m# = m# - twopi’ in the (-pi, pi) interval 

’Compute the initial approximation: 
g# = 1 / (4 * ec# ♦ .5) 
a# = (1 - ec#) * g# 
b# = .5 * m# * g# 

arg# = ABS(b#) + SQR(b# * b# ♦ a# * a# * a#) 
z# = arg# ‘ third 
IF b# < 0 THEN z# = -z# 
s# = z# - a# / z# 

s# = s# - .078 * s# " 5 / (1 ♦ ec#) 

ea# = m# ♦ ec# * s# * (3 - 4 * s# * s#) ’ abs(dea/ea) <= 0.002 
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'Refine the approximation: 



f 2# 


= 


ec# * 


SIN(ea#) 


f 3# 


= 


ec# * 


C0S(ea#) 


fO# 


= 


ea# - 


f 2# - m# 


f 1# 


= 


1# - 


f 3# 


f 4# 


= 


-f 2# 




f 5# 


= 


-f 3# 





'This next correction term is sufficient if ecc < 0.25 
d# = -fO# / fit 

'Use more of the following terms to increase accuracy 
d# = -fO# / (fl# ♦ .5# * f 2# * d#) 

df = -fO# / (fit ♦ d# * (.5# * f 2# ♦ d# * r3 * f3#)) 

d# = -fO# / (fl# + d# * (.5# * f 2# ♦ d# * (r3 * f3# ♦ d# * r4 * f4#))) 

d# = -fO# / (fl# ♦ d# * (.5# * f 2# + d# * (r3 * f3# ♦ d# 

» (r4 * f 4# + d# * r5 * f5#)))) 
ea# = ea# + d# 

END IF 
END SUB 



>*************************************************** *x********x*»*«xx*xx**xxx* 

SUB output4 (ky., yr&, moft, dafc, time. of. day, It, In, ht , el, az , rng , _ 
langle, head, 8ta.id$, sat.idS) 

CONST pi = 3.141592653589793#, rad = pi / 180#, tp = pi ♦ pi 
c$ = C HRS (34) 



SELECT CASE X 1 /. 
CASE 0 



PRINT 


n 


height 


PRINT 


" look" 




PRINT 


"year mo da hr mn sec lat 


long (km) 


PRINT 


"azim range angle head" 




PRINT 

CASE 1 

CALL hr. min. 


sec(time . of .day# , hry., miny,, 


sec#, 1) 



It = It / rad 

In = dmod(ln / rad, 360#) 



elev 



la = 
hd = 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
CASE ELSE 
PRINT 
PRINT 
PRINT 

END SELECT 

END SUB 



langle 
head / rad 

USING "#### ## ## ## ## ##.# yrk ; mok ; daA; hr*/.; mm*/.; 
USING "###.## ###.## #####.# ###.## It; In; ht ; el; 
USING "###.## ##### ##.## ###.##"; az; rng; la; hd 
#3, USING "!*! \k\ M ; c$; sta.idS; c$; c$; sat.idS; c$ ; 

#3, USING M !####! !##! M ; c$; yrk; c$; cS; mok ; c$ ; 

#3, USING “!##! !##! "; c$ ; da&; c$ ; c$; hry.; c$ ; 

#3, USING "!##! !##.#! "; c$; miny.; c$; c$; sec; c$ ; 

#3, USING "!###.##! !###.##! " ; cS; It; c$ ; c$; In; c$ ; 

#3, USING "!#####.#! !###.##! c$ ; ht ; c$; c$; el; c$; 
#3, USING "!###.##! !#####! " ; c$; az ; c$; c$; rng; c$ ; 

#3, USING "!##.##! !###.##!"; c$; la; c$; c$; hd ; c$ 



"Error in output4 subroutine. Must have 0 <= ky. <= 1." 

"ky. = ky. 



M . 



sec ; 



SUB pos. update (ll , elem#(), time#, sat.x#(), sat.xd#(), elemu#(), r# , ea#) _ 
STATIC 

> 09-07-88. Rev. 03-09-89. 

> Case 1: Initialize orbital elements for time of epoch. 

> Case 2: Update orbital elements and compute in-plane position k velocity. 

» Simplest of NAVSPASUR methods. Accurate to 10-15 n.mi. within 20 days 
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of epoch. 



' elem(l) 
' elem(2) 
' elem(3) 
' elem(4) 
' elem(5) 
' elem(6) 
' elem(7) 
' elem(8) 
' elem(9) 



amajor 

ecc 

tzero 

aincl 

anode 

peri 

motion 

manomaly 

decay 



semi-major axis 
eccentricity 

time of perifocal passage 

inclination (radians) 

longitude of ascending node (radians) 

argument of perifocus (radians) 

mean motion (revolutions/day) 

mean anomaly (radians) 

1st decay constant (radians/minute~2) 



CONST pi = 3.141592653589793#, rad = pi / 180#, twopi = pi ♦ pi 
CONST c = .00081197385# ' 0.75*j2 

CONST twothird = 2# / 3# 

CONST ae = 6378.14# ' earth's equatorial radius (km.) 

CONST GE = 398600.5# ' geo. grav. const, (km) ~3/(sec) “2 

k = 60# * SQR(GE / ae) / ae ' (eru)- (3/2)/min 
CONST we = 1.002737909350795# * twopi / 1440# ' sidereal period in radians/min 



SELECT CASE lX 
CASE 1 

elem#(l) = (k / elem#(7)) ~ twothird 

a2# = elem#( l) * 2 

ecc2# = 1# - elem#(2) “ 2 

ci# = C0S( elem#(4) ) 

c2i# = ci# * ci# 

elem#( 1 ) = elem#(l) * (1# ♦ c * (3# * c2i# - 1#) / _ 
(a2# * ecc2# * 1.5#)) ~ twothird 
elem#(3) = time# - elem#(8) / elem#(7) 
aecc22# = (elem#(l) * ecc2#) “ 2 
peridot# = c * (5# * c2i# - 1#) / aecc22# 
nodedot# = -2# * c * ci# / aecc22# 
adot# = -2# * twothird * elem#(l) * m2# / elem#(7) 

FOR it = 1 TO 9 

elemu#(iX) = elern#(iX) 

NEXT 
CASE 2 

ma.mO# = (elem#(9) * time# + elem#(7)) * time# 
ma# = ma.mO# + elem#(8) 

elemu#(6) = elem#(6) ♦ peridot# * ma.mO# 

elemu#(5) = elem#(5) ♦ nodedot# * ma.mO# - we * time# 

elemu#(l) = elem#(l) + adot# * time# 

CALL kepler(ma#, elem#(2) , eat) 
sO# = SIN(ea#) 
cO# = C0S(ea#) 
elemu#(8) = ma# 

elerau#(3) = time# - ma# / elem#(7) 



' in-plane position and velocity 
sat.x#(l) = elemu#(l) * (cO# - elem#(2)) 
sat.x#(2) = elemu#(l) * SQR(ecc2#) * sO# 
sat.x#(3) = 0# 

r# = elemu#(l) * (1# - elem#(2) * cO #) 

ed# = elemu#(l) * elem#(7) / r# 

sat.xd#(l) = -elemu#(l) * ed# * sO# 

sat.xd#(2) = elemu#(l) * ed# * SQR(ecc2#) * cO# 

sat.xd#(3) = 0# 

' equatorial rectangular coordinates 

CALL euler (elemu(6) , elemu(5), elem(4) , sat.xO, sat.xdO) 

' correct velocity for earth's rotation 
sat.xd#(l) = 8at.xd#(l) + we * sat.x(2) 
sat.xd#(2) = sat.xd#(2) - we * sat.x(l) 
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CASE ELSE 
PRINT 

PRINT “Error in pos . update. First argument must be 1 or 2 only." 
PRINT "Argument = 1% 

PRINT 

END SELECT 
END SUB 



J ***************************************************** ******* *************** 
SUB rise. set (zv#(), sta.R#(), wet, tO#, eat, elemu#()) 

1 Compute eccentric anomaly for the time of rise or set of a satellite 
' 02-04-89. Rev. 02-08-89 C 1700. 

* NOTE: To use this routine, the time of culmination MUST be computed first, 

* The entering ea# argument should be ea.cul# -/♦ 0.2 radians (about 11 
' degrees) for rise/set. 

} See the introductory comments in SUB culminate for nomenclature and 
' notation. 

CONST pi = 3.141592653589793#, rad = pi / 180#, twopi = pi ♦ pi 

amaj# = elemu#(l) 
ecc# = elemu#(2) 
capt# = elemu#(3) 
incl# = elemu#(4) 
node# = elemu#(5) 
peril = elerau#(6) 
mot# = elemu#(7) 

DIM pv#(3), pvd#(3) , qv#(3), qvd#(3) , rvecd#(3) , rho#(3) 

cosw# = C0S(peri#) 
sinw# = SIN(peri#) 
cosi# = C0S(incl#) 
sini# = SIN(incl#) 



corr# = 1 
DO 

cosea# = C0S(ea#) 
ecosea# = ecc# * cosea# 
sinea# = SIN(ea#) 
esinea# = ecc# * sinea# 



gamma# = node# - we# * ((ea# - esinea#) / mot# ♦ capt# - t0#) 
gammad# = -we# * (1# - ecosea#) / mot# 



cosg# = C0S(gamma#) 
sing# = SIN(gamma#) 



rtecc# = amaj# * SQR(l# - ecc# * ecc#) 



xw# = amaj# * (cosea# - ecc#) 
xwd# = -amaj# * sinea# 
yw# = rtecc# * sinea# 
ywd# = rtecc# * cosea# 



pv#( 1) 
pv#(2) 
pv# (3) 
qv#( 1 ) 
qv#( 2) 
qv#( 3) 



cos w# * cosg# - sinw# * sing# * cosi# 
cosw# * sing# + sinw# * co sg# * cosi# 
sinw# * sini# 

-sinw# * cosg# - cosw# * sing# * cosi# 
-sinw# * sing# ♦ cosw# * cosg# * cosi# 
cosw# * sini# 
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pvd#(l) = -pv#(2) * gammadf 
pvd#(2) = pv#(l) * gammadf 
pvdf (3) = Of 

qvd# ( 1 ) = -qv#(2) * gammadf 
qvd#(2) = qvf(l) * gammadf 
qvdf (3) = Of 



FOR iX = 1 TO 3 

rhof(iX) = xw# * pvf(iX) ♦ yw# * qvf(iX) - sta.Rf(iX) 
rvecdf(iX) = xwd# ★ pvf(iX) ♦ xw# * pvdf(iX) ♦ ywd# * 
+ yw# * qvdf(iX) 



NEXT 



qvf(iX) . 



fe# = dotf(zv#(), rhof()) 
fed# = dot#(zv#(), rvecdfO) 

corrf = fef / fed# 
eaf = eaf - corrf 



LOOP UNTIL ABS( corrf ) < .OOOOlf 
END SUB 



>*************************************************** c****************»c********** 

SUB sat. head (rf, sat.x#(), sat.xd#(), headf) 1 01-18-89. Rev. 03-12-89. 

'satellite heading 

CONST pi = 3.141592653589793#, twopi = pi + pi 

shead# = (sat.x#(l) * sat.xd#(2) - sat.xf(2) * sat.xdf(l)) / r# 
chead# = sat.xd#(3) - (sat.x#(3) * dot# (sat . xf ( ) , sat.xdfQ)) / r# * 2 
headf = atan2# ( shead# , cheadf) 

IF headf < 0# THEN head# = headf ♦ twopi 

END SUB 



'ft***************************************************************************** 

SUB sat . sphere . coord (sat.x(), r, t j , tm, sat.dc, sat. In, sat.ra) 

'compute spherical coordinates 



CONST pi = 3. 141592653589793f , rad = pi / 180f, twopi = pi ♦ pi 



sat.dc = asin(sat . x(3) / r) 
sat. In = atan2(sat . x(2) , sat.x(l)) 
ru = 15# * gmst(tj, tm) * rad 
sat.ra = dmod(sat.ln + ru, twopi) 



1 satellite declination 
' satellite longitude 
' Greenwich mean sidereal time 
' satellite right ascension 



END SUB 



>************************************************ ******** ********* ************* 
SUB sta.p. coord (ltf, In#, htf, x#Q) ' 08-06-88. Rev 02-04-89. 

' convert station, lat , long k height (ft.) from geographic 
' to geocentric rectangular - x(.) 



CONST nmi.ft = .3048# / 1852# 
CONST km. ft = .3048# / 1000# 
CONST ae = 6378.14# 

CONST fl = If / 298.257# 

CONST e.sqr = (2# - fl) * fl 



' n. mi. /foot 
9 km. /foot 

' earth's equatorial radius (km.) 
' earth's flattening factor 
' earth's ellipticity squared 



ht.eruf = htf * km. ft / ae 

sit# = SIN(ltf) 

clt# = COS(ltf) 

sin# = SIN(lnf) 

cln# = C0S(ln#) 

nl s 1# / SQR(lf - e.sqr * sit * sit) 
gif = nf + ht . eruf 
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g2f = nf * (if - e.sqr) ♦ ht erul 



Radius vector for flattened earth 
xf (1, 1) = gl# * clt# * clnf 

xf(2, 1) = gl# * cl tf * elnf 

xf (3 , 1) = g2f * sltf 

Spherical earth unit vectors: 

1 1) local perpendicular vector 
xf (1, 2) = cltf * clnf 
xf (2 , 2) = cltf * s Inf 
xf ( 3 , 2) = sltf 

’2) east vector 
xf(l, 3) = -8 Inf 

xf (2 , 3) = clnf 

xf (3 , 3) = Of 

*3) north vector 
xf (l , 4) = -sltf * clnf 
xf(2, 4) = -sltf * slnf 

xf(3, 4) = cltf 

END SUB 

J*****XX*****XX*x******X********************x****x*******X******xxxxxxxxx»««x* 

SUB sta.sat (sta.xfQ, sat.xf(), rngf , azf, elf, langlef) 

* 08-05-88. Rev. 08-31-88. 

'compute station-satellite relations. Range, azimuth, elevation and satellite 
' "look-angle * 1 (angle between satellite Bub point and station). 

CONST pi = 3 . 14l592653589793f , rad = pi / 180f, twopi = pi ♦ pi 

'rho() = sat coordB minus sta coords, rhodots are the rho components projected 
' in the station local north-east-radial system, 

DIM rhof (3) , rhodotsf(3) 

rngf = Of 
FOR iX = 1 TO 3 

rf = sat.xf(iX) - sta.xf(iX, l) 

rhof(iX) = rf 

rngf = rngf ♦ rf * rf 

NEXT 

rngf = SQR(rngf) 

FOR iX = 1 TO 3 

rhof(iX) = rhof(iX) / rngf 

NEXT 

FOR iX = 1 TO 3 
af = Of 

FOR jX = 1 TO 3 

af = af ♦ rhof(jX) * sta.xf(jX, iX ♦ l) 

NEXT 

rhodotsf(iX) = af 

NEXT 

azf = atan2f ( rhodot Bf (2) , rhodotBf(3)) / rad 
IF azf < Of THEN azf = azf ♦ 360f 

af s rhodot8f(l) 

elf = ATN(af / SQR(lf - af * af)) / rad 

af = Of 
bf = Of 

FOR iX = 1 TO 3 
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c# = sat.x#(iX) 

a# = a# + rho#(iX) * c# 

bt = bf + ct * c# 

NEXT 

langle# = acos#(a# / SQR(b#)) / rad 
END SUB 

>***************************************************************** ************* 
SUB time. update (kX, delta. tf, tm#, jd#, jd#(), tj#, time. of. day, mJt, dfc, _ 
yt t hf, mo$, day$) 

CONST twentyfour = 24#, sixty = 60#, jul.cent = 36525# 

t# = delta. t# 'delta. t# in hours if kX = 0, t# in hours. 

IF kX > 0 THEN t# = t# / sixty 'delta. t# in minutes if kX > 0, t# in hours, 
tm# = tm# ♦ t# 

jd#(2) = jd#(2) ♦ t# / twentyfour 
jd# = jd#(l) 4 jd#( 2) 

tj# = tj# 4 t# / twentyfour / jul.cent 
time. of. day = time. of. day 4 t# 

IF time. of. day >= twentyfour THEN 

time. of. day = time. of .day - twentyfour 
jd#(l) = jd#( 1) 4 l#; j d# ( 2 ) = j d# C 2 ) - 1# 

CALL jd .to .date( jd#( ) , m&, d&, y&, h# , mo$ , day$) 

END IF 

IF time. of .day < 0# THEN 

time. of. day = time. of. day 4 twentyfour 
j d # ( 1 ) = j d# ( 1 ) - 1#: j d# ( 2 ) = jd#(2) + 1# 

CALL jd . to .date( jd#() , m&, d&, yk, h# , mo$ , day$) 

END IF 

END SUB 

SUB yrmoday (v$, yr&, mot, dy&) ' 03-21-88. Rev 03-21-88. 

' Convert yyyy.mmdd to yrfc=yyyy, mot=mm, and dyfc=dd 

IF INSTR(v$, ".") <> 5 OR LEN(v$) <> 9 THEN 
PRINT 

PRINT "Error in yrmoday. v$ must be of the form yyyy.mmdd" 

PRINT "v$ = "; v$ 

PRINT 

ELSE 

yrk = VAL (LEFT$ ( v$ , 4)) 
mot = VAL(MID$(v$, 6, 2)) 
dyk = VAL(RIGHT$(v$, 2)) 

END IF 
END SUB 
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APPENDIX E. THE HIGH ACCURACY KINEMATICS ROUTINE. 



This appendix contains a listing of the alternate pos. update routine. It is a BASIC ver- 
sion of the highly accurate PPT2 subroutine provided by NAVSPASUR [Ref. 1] with their 
SHOWALL program. It is believed that it has been accurately translated, however any errors 
are solely those of the author. Known discrepacies are the difference in the WGS-72 con- 
stants used in SHOWALL and the IAU 1976 constants used here [see section III. SOURCES, 
in this report for details]. 

*********** ************** *********************************************** 

SUB pos. update (indX, elem#(), timet, sat.x#(), sat.xd#(), elemu#(), _ 
rt, eat) STATIC Rev. 09-26-89 0 0900 

' Case 1: Initialize orbital elements for time of epoch. 

9 Case 2: Update orbital elements and compute in-plane position t velocity. 

' Comments from original FORTRAN version: 

9 implicit real*8(a-h , o-z) 

>*** model r ef erence-Astronomical Journal 64, No 1274 , Nov . , * 59 
>*** solution of the problem of artificial satellite theory 
>*** with out drag... Dirk Brouwer. 

9 840329 confirmed use of w=u x v vice w=u x vel 

9 840329 inserted cdh,sdh, cdt,sdt re osc vice mean els at t 

* for use only with dcext2 , del ,dcia 

9 frame for u,v,w is the rotated inertial w/o coriolis 

9 830127 new partials re ( a , x , y , z , p , q) 

* 830127 improved (polynomial) model for singularity in Brouwer 
DIM a( 25 ) 

DIM f (25 ) , osc(10), kf(10), cf(10), bs(3, 4), u(3), v(3), w(3), vel(3) 

DIM ar (6 , 8), br(3, 6), cr(3, 6), dv(3) 



**** The commented beta value is the WGS-72 value , however using this 
>*** value, one obtains a value for a(9) different then what we had 
>*** been using, thus causing a problem with time . therefore to retain 
>*** previous value of a(9),one solves for beta using the old value 
'*** of a(9) . 

**** betal=398600 . 5d0 

> betal = 398597.62579588#, erkm = 6378.135#, flat = 298.26# 

1 k = .0743669161# 



' IAU 1976: 

CONST twopi = 6.283185307179586# 

CONST we = 1.002737909350795# * twopi 
CONST ae = 6378.14#, erkm = ae 
CONST betal = 398600.5# 

k = 60# * SQR(betal / ae) / ae 
CONST flat = 1# / 298.257# 



/ 1440# ' sidereal period in radians/min 
1 earth's equatorial radius (km.) 

9 geo. grav. const, (km) "3/(360) “2 
1 (eru)~(3/2)/min 
• earth's flattening factor 



' k-terms [Note: These may or may not agree with IAU 1976] 

CONST c20 = -.0004841605#, c30 = .00000095958# 

CONST c40 = .00000055199#, c50 = .000000065875# 



CONST twotrd = .666666666667#, fortrd = 1.333333333333# 
CONST tentrd = 3.333333333333#, ar = 48 * 0# 



'IF ind^ <> 1 THEN GOTO secondcase 
SELECT CASE indX 
CASE 1 



kzX = 1 

a( 1) = 0#' [Note: = 85-01-01 in SHOWALL. Set to zero here.] 
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a( 3) = -.5# * c20 * SQR(5i) 
a(4) = c30 * SQR(7f) 
a(5) = .375 * c40 * SQR(9f) 
a( 6) = c50 * SQR(llf) 

* two pi 
a(7) = twopi 

a( 16) = It / a(7) 

* hergs/day , secs/herg, mins/herg 
a(9) = erkm * SQR(erkm / betal) 

a(8) = 86400# / a(9) 
a( 17) = 1440# / a(8) 

' we(rad/day), we- 2 pi, we(rad/herg) 

a( 1 1 ) = 6.3003880987# 
a( 10) = a( 11) - a(7) 
a( 12) = a(l 1 ) / a(8 ) 

, earth flattening 

a( 13) = (2# - 1# / flat) / flat 

* sm/er, km/er, nm/er 
a(20) = erkm / 1.609344# 
a(21) = erkm 

a(22) = erkm / 1.852# 

* deg/rad 
a(23) = 360# / twopi 

* range rate/er/herg to cycles/second - conversion 
a(24) = a(21 ) * 216980000# / (a(9) * 299792.5#) 

* fence plane displac from earth center 
a( 25) = .0031# 



Patch to interface PPT2 and SATSTA variables : 



f (1) 
f ( 2) 
f ( 3) 
f (4) 
f ( 5) 
f ( 6) 
f (7) 
f(8) 
f ( 9) 



elera(8 ) 
elem(7) 
elem(9) 

0 # 

elem(2) 

elem(6) 

elem(5) 

C0S( elem(4) ) 



'mean anomaly 
a(l7) 'mean motion 
a(l7) * 2'first decay constant 
Second decay constant 
' eccentricity 
'arg of perigee 
'node 

'cosine inclination 



= 0 # 



epoch in PPT2, time from epoch here (? but check ?) 



' Used only once to “recover" the remaining elements 
hi = 0# 

'*** f (8) is cosine of inclination 
t9 = f (8) * f (8) 
tl = 1# - 5# * t9 

' t2=( 1 . 0d0-exp(-100 . 0d0*t 1 ~2) )/t 1 

' revised handling of divisor per navspasur o2t memo nov , 83 

beta = 100# / 2# “ 11 

p3 = beta * tl ~ 2 

pi = EXP (-p3) 

p2 = 1# ♦ pi 

plex = pi * 2 

p4 = 1# 

p3ex = - . 5# * p3 
FOR nX = 2 TO 13 

IF nX <= 11 THEN p2 = p2 * (1# ♦ plex) 

plex = plex * plex 

p4 = p4 ♦ p3ex 

p3ex = -p3 * p3ex / (nX ♦ 1) 

NEXT 

t2 = p2 * p4 * beta * tl 
cio2 = SQR( . 5# ♦ .5# * f(8)) 
sio2 = SQR(.5# - .5# * f (8) ) 
tl2 = 3# * t9 - 1# 

tl3 = t9 * t9 

t 14 = t2 * t 13 

t 15 = 1# - t9 
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'*** f(5) is eccentricity 

esqO = 1(5) * 1( 5) 

eta20 = 1# - esqO 

etaO = SQR(eta20) 

t6 = etaO * eta20 

t7 = eta20 * eta20 

t8 = etaO 4 if / (if + etaO) 

1(13) = 1(2) * (-tvotrd) 
elem(l) = f(l3) 

>+*+ 1(2) is basically mean motion but also contains the 

'*** secular term in mean anomaly, therelore 1(13) is (almost) the 

'*** semi major axis. 

This loop is to remove the contribution ol the additional 
>+*+ term in mean motion lor improving values ol semi-major axis 
'*** and seculars. 

' *** a(3) to a(6) are 'k* terms wgs 72 

FOR i% = 1 TO 5 

1 s 13 = 1(13) * 1(13) * t7 

Isll = a(3) / lsl3 

Is 13 = a(5) / (1 s 13 * Is 13) 

tt = If 4 1 . 5f ♦ Isll * etaO * tl2 

uu = -I5f + 16# * etaO 4 25f * eta20 

uu = uu + (30# - 96# * etaO - 90# * eta20) * t9 

uu = uu + (105# + 144# * etaO ♦ 25# * eta20) * tl3 

tt = tt + .09375# * Isll 2 * etaO * uu 

tt = tt ♦ .9375# * lsl3 * etaO * esqO * (3# - 30# * t9 + 35# * t!3) 

1(13) = (tt / 1(2)) * twotrd 

NEXT 

>*** f(ll) will be rate ol change ol argument ol perigee 

uu = -35# + 24# * etaO + 25# * eta20 

uu = uu + (90# - 192# * etaO - 126# * eta20) * t9 

uu = uu + (385# + 360# * etaO + 45# * eta20) * tl3 

tt = -1.5# 4 Isll * tl ♦ .09375# * Isll “ 2 4 uu 

uu = 21# - 9# * eta20 + (-270# + 126# * eta20) * t9 

uu = uu 4 (385# - 189# * eta20) * tl3 

1(11) = tt + .3125# * Is 13 * uu 

>*** 1(12) will be rate ol change ol right ascention 

uu = -5# + 12# * etaO + 9# * eta20 

uu = uu - (35# + 36# * etaO 4 5# * eta20) * t9 

tt = -3# * Isll * 1(8) 4 .375# * Isll * 2 * 1(8) * uu 

1(12) = tt 4 1.25# * Is 13 4 1(8) * (5# - 3# * eta20) * (3# - 7# * t9) 

’*** l(l5)=sine inclination ,f ( 14)i8 a dot,l(16) is eccen.dot 
1(15) = SQR(tlS) 

1(14) = -lortrd * 1(3) * 1(13) / 1(2) 

1(16) = 1(14) * f (5) ♦ et a20 / 1(13) 

ql = 1(2) 4 1(13) * 1.5# 

1(11) = 1(11) / ql 

1(12) = 1(12) / ql 

Is 12 = a(4) / (a(3) 4 f(l3) 4 e ta20) 

1 s 13 = 1 s 13 / Isll 4 tentrd 

Is 14 = a(6) / (a(3) ♦ f (13) * 3 4 e ta20 4 t 7) 

ql = .125# 4 (fall 4 (l# - 11# * t9 - 40# 4 tl4) - 1 sl3 . 

4 (1# - 3# 4 t9 - 8# 4 t 14 ) ) 

p5 = 1# 4 t2 4 (8# 4 t9 4 20# 4 t 14) 

p2 = 1# + 2# 4 p 5 

q2 = .125# 4 esqO 4 f(8) 4 (f 8 ll 4 (l# + io# ♦ p 5) - f s 13 ♦ p2) 

p2 = .46875# 4 p 2 * 1(5) 4 1(8) 4 1(15) 4 (4# 4 3# 4 esqO) * fsl4 

p3 = Is 14 4 (1# - 9# 4 t 9 - 24# * 1 14) 
q5 = .25# 4 (1 8 12 + .3125# 4 ( 4 # 4 3 # ♦ esqO) 4 p 3) 
p3 = .15625# 4 1(5) 4 f(i5) 4 p 3 

p4 = .030381944# 4 f(5) 4 ial4 * ( 1 # - 5# * t9 - 16# 4 t 14) 

p5 = .060763889# 4 f(5) 4 e8 q0 * 1(8) ♦ 1(15) 4 Ial4 4 (l# 4 4 # * p5) 

vl el = 1(5) 4 e ta20 ♦ ql 

vlhli = -f CIS) 4 q2 

tt = ( t6 - 1#) * ql - q 2 4 25# 4 e sqO * t!4 4 t 9 ♦ t 2 4 (fsll - . 2 # * fsl3) 
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vlsl = tt 



- .0625# * esqO * (fall * (It - 33# * 

- Is 13 * (1# - 9# * t9 - 40# * t 14) ) 
vle2 = eta20 * 1(15) * q5 

vlh2i = 1(5) * f (8) * q5 ♦ f(l5) * p2 

vls2 = 1(5) * 1(15) * (t8 ♦ 1(8) / (1# ♦ 1(8))) 

♦ (11# ♦ 3# * esqO - 3# * t6) * p3 ♦ 
vl©3 = -3# * 1(5) * ©ta20 * 1(15) * p4 
vlh3i = -esqO * 1(8) * p4 - 1(15) * p5 
v!b3 = 1(15) * (3# * t6 - 3# - 2# * ©sqO - ©sqO 
- (1# - 1(8)) * P 5 

vll2 = vl©2 ♦ 3# * f (5) * eta20 * p3 

} precomps lor the partials 

'*** 1(6) is argument ol p©rig©0,l(7) is right 

sinhO = SIN(1(7)) 

coshO = C0S(1 (7) ) 

sintO = SIN(1(6) ♦ 1(7)) 

costO = C0S(1(6) ♦ 1(7)) 

’ kepler 1 s law 

dadm = -twotrd * 1(13) / 1(2) 

1 variations thru the seculars 

secul = a(3) / (1(13) * eta20) * 2 

dgddm = 3# * secul * tl / 1(13) * dadm 

dgdde = -6# * secul * tl * 1(5) / eta20 

dgddi = -15# * secul * 1(8) * 1(15) 

dhddra = 6# * secul * 1(8) / f(l3) * dadm 

dhdde = -12# * secul * f(8) * 1(5) / eta20 

dhddi = 3# * secul * 1(15) 

dtddx = (dgdde ♦ dhdde) * costO 

dtddy = (dgdde ♦ dhdde) * sintO 

dtddp = (dgddi ♦ dhddi) * 2# * coshO / cio2 

dtddq = (dgddi + dhddi) * 2# * sinhO / cio2 

dhddx = dhdde * costO 

dhddy = dhdde * sintO 

dhddp = dhddi * 2# * coshO / cio2 

dhddq = dhddi * 2# * sinhO / cio2 

' variations thru m2 

daddml = lortrd * 1(3) * (1(13) / 1(2) - dadm) / 
daddm2 = -lortrd * 1(13) / f(2) 

dedde = -fortrd * (l# - 3# * esqO) * 1(3) / 1(2) 
deddx = dedde * costO 
deddy = dedde * sintO 

deddml = lortrd * f(5) * eta20 * 1(3) / 1(2) “ 2 
deddm2 = -fortrd * f(5) * eta20 / f(2) 

FOR iX = 1 TO 9 

elemut(iX) = elem#(iX) 

NEXT 

* secondcase : 

CASE ELSE 

tm = time# / a(17) 'Convert minutes to Hergs 
hi = tm - 1(9) 

>*** compute position as a lur.ction ol time (in 
osc (9) * hi * (1(2) ♦ hi * (1(3) ♦hi * 1(4))) 
osc(3) = dmod(l(l) ♦ osc(9), twopi) 
tt = 1(5) ♦ 1(16) * hi 
IF tt < 0# THEN tt = 0# 

IF tt > .99999# THEN tt = .99999# 

08c(4) =• tt 
esq = 08c(4) ~ 2 
eta2 = 1# - esq 
eta = SQR(eta2) 

CALL kepler(osc(3) , osc(4), c3) 
osc(l) = C0S(c3) 
cl = 1# - osc(4) * osc(l) 
osc(l) = (osc(l) - 08C(4)) / cl 



t9 - 200# * t 14) _ 

* q5 _ 

( 1 # - 1 ( 8 )) * P 2 

* 1(8) / (1# ♦ 1(8))) * p4 _ 

ascension 



1 ( 2 ) 



call line) 
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os c(2) = eta * SIN(c3) / cl 

ts4 = If ♦ osc(4) * osc(l) 
cf 3 = 0# 

IF kz% <> 0 THEN cf 3 = a(l2) 
tt = f ( 13) ♦ f ( 14) * hi 

IF tt < It THEN tt = It 

08C(8) = tt 
osc(7) = f(8) 
fsi = f ( 15) 

osc(6) = dmod(f(7) ♦ f(12) * osc(9), twopi) 

os c(5) = dmod (f (6) ♦ f(ll) * osc(9), twopi) 

r = osc(8) * eta2 / ts4 

osc(6) = osc( 6) - cf 3 * (tm - a(l)) 

wl = COS (oscC 5) ) 

w2 = SIN(osc(5)) 

w7 = C0S(osc(6)) 

w8 = SINCosc ( 6) ) 

agda = Ot 

agde = Ot 

agdi = Ot 

agdl = Ot 

agdg = Ot 

agdh = Ot 



w3 


= wl 


♦ 


wl 


- w2 


♦ 


w2 


w4 


= 2t 


* 


wl 


* w2 






w5 


= w3 


* 


wl 


- w4 


* 


w2 


w6 


= w4 


* 


wl 


♦ w3 


♦ 


w2 



w9 = osc(l) * osc(l) - osc(2) * osc(2) 

w 10 = 2t * osc(l) * osc (2) 

wll = w3 * osc(l) - w4 * osc(2) 

wl2 = w4 * osc(l) ♦ w3 * osc(2) 

wl3 = w3 * w9 - w4 * wlO 

wl4 = w4 * w9 ♦ w3 * wlO 

w!5 = wl3 * osc(l) - wl4 * osc(2) 

wl6 = wl4 * o 8 c ( 1 ) ♦ wl3 * osc (2) 

wl7 = atan2(osc(2) , osc(l)) ♦ osc(4) * osc(2) - dmod(osc(3), twopi) 
w 1 8 = C0S(osc ( 3) ) 
w 19 = SIN(osc (3) ) 

w20 = osc(l) * (3t + osc(4) * osc(l) * (3t ♦ osc(4) * osc(l))) 
w21 = 3t * w 14 ♦ 3t * osc (4 ) * wl2 ♦ osc(4) * wl6 
w22 = ts4 * (it ♦ ts4) / eta20 

osc(8) = osc(8) * (it ♦ fsll / eta20 * (tl2 * (ts4 3 - t6) _ 

♦ 3t * tl5 * ts4 3 * wl3)) ♦ agda 
de = vlel * w3 ♦ vle2 * u2 ♦ vle3 * w6 

di = sio2 ♦ . 5t * cio2 * ( . 5t * fsll * f (8) * f(15) * (3t * wl3 ♦ 3t * f(5) _ 

♦ wll + f (5) * wl5) - de * f (5) * f (8) / f ( 15 ) / eta20) ♦ .5t * cio2 * agdl 
de = osc(4) + . 5t * fsll * ( 1 12 * (w20 ♦ f(5) * t8) ♦ 3t * tl5 * (w20 ♦ f(5)) . 

♦ wl3 - eta20 * tl5 * (3t * wll + wl5)) ♦ de ♦ agde 

dh = . 5t / cio2 * ( - . 5t * fsll * f(8) * f(15) * (6t * wl7 - w21) ♦ vlhli * w4 _ 

♦ vlh2i * wl ♦ vlh3i * w5) ♦ .5t * agdh / cio2 
dl = - . 25t * t6 * fsll * (2t * 1 12 ♦ ( w22 ♦ It) * osc(2) ♦ 3t * tl5 * ((it _ 

- w22) * wl2 ♦ ( . 333333333t + w22) * wl6)) 
os = osc(3) + osc(5) ♦ osc(6) - dl * f ( 5 ) * (t8 - It) / t6 - .25t * fsll * (6t _ 

♦ wl7 * (tl ♦ 2t * f (8) ) - w21 * ( t 1 ♦ 2t 4 2t * f(8))) _ 

♦ vlsl * w4 ♦ vls2 * wl + vls3 * w5 ♦ agdg 

dl = dl ♦ etaO * (vlel * w4 - vll2 * wl - vle3 * w5) ♦ agdl 
esq = de * de ♦ dl * dl 
osc(4) = SQR(esq) 

os c (3) = atan2(de * wl9 ♦ dl * wl8, de * wl8 - dl * wl9) 
osc(7) = It - 2t * (di * di ♦ dh * dh) 
fsi = SQR( It - osc( 7) * 08c(7) ) 

os c(6) = atan2(di * w8 ♦ dh * w7, di * w7 - dh * w8) 
os c(5) = os - osc(3) - osc(6) 
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eta2 = It - esq 
eta = SQR(eta2) 

CALL kepler(osc(3) t osc(4), c3) 
osc(l) = C0S(c3) 
cl = It - osc (4) * osc(l) 
osc(l) = (osc(l) - osc(4)) / cl 
osc(2) = eta * SIN(c3) / cl 

ts4 = 1# ♦ osc(4) * osc(l) 
r = osc(8) * eta2 / ts4 
w7 = COS (osc(6) ) 

¥8 = SIN(osc( 6) ) 

Ml = C0S(osc(5) ) 

w2 = SIN(osc(5) ) 

ubs = w2 * osc(l) ♦ wl * osc(2) 

ubc = wl * osc(l) - w2 * osc(2) 

u(l) = m! * ubc - w8 * ubs * osc(7) 

u(2) = w8 * ubc ♦ m! * ubs * osc(7) 

u(3) = ubs * fsi 

v(l) = -w7 * ubs - w8 * ubc * osc(7) 

v(2) = -»8 * ubs ♦ w7 * ubc * osc(7) 

v( 3) = ubc * fsi 

w(l) = w8 * fsi 

w(2) = -w7 * fsi 

w(3) = osc(7) 

ts2 = eta * SQR(osc(8)) 

FOR iX = 1 TO 3 

vel(iy.) = (osc(4) * osc(2) * u(iX) + ts4 * v(iX)) / ts2 

NEXT 

vel(l) = vel(l) ♦ cf3 * r * u(2) 
vel(2) = vel(2) - cf3 * r * u(l) 

' Patch for PPT2 to SATSTA output 
elemu(l) = osc(8) * semimajor axis 
elemu(2) = osc (4) 1 eccentricity 

’elemu^) 2 'time of last perigee passage 

elemu(4) = acos (osc (7) )' inclination 

elemu(5) = osc(6)'node 

elemu(6) = osc(5) , arg of perigee 

elemu(8) = osc(3)'mean anomaly 

mat = elemu(8) 

CALL kepler(mat, elem#(2) , eat) 
sOf = SIN(eaf) 
cOf = COS(eaf) 

1 elemut(8) = mat 

elemu#(3) = timet - mat / elemt(7) 

' in-plane position and velocity 
ecc2t = It - elemut (2) * 2 
sat.xt(l) = elemut(l) * (cOt - elemt(2)) 
sat.xt(2) = elemut(l) * SQR(ecc2t) * sOt 
sat.xt(3) = Ot 

rt s elemut(l) * (It - elemt(2) * cOt) 

edt = elemut(l) * elemt(7) / rt 

sat.xdt(l) = -elemut(l) * edt * sOt 

sat.xdt(2) * elemut(l) * edt * SQR(ecc2t) * cOt 

sat . xdt(3) = Ot 

' equatorial rectangular coordinates 

CALL euler(elemu(6) , elemu(5), elem(4) , sat.x(), sat.xdQ) 

END SELECT 
END SUB 
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